Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


ATIQN  PAGE 


lb.  RESTRICTIVE  MARKINGS 


4.  PERFORMING  ORGANIZATION  REPORT  N8jM|E| 


6a.  NAME  OF  PERFORMING  ORGANIZATION 
University  of  Oslo 


6c.  ADDRESS  {City,  State,  and  ZIP  Code) 


FFICE  SYMBOL 
(If  applicable) 


Faculty  of  Science,  PO  Box  1032  Blindern 
N-0316  OSLO,  NORWAY 


8a.  NAME  OF  FUNDING  /  SPONSORING 
ORGANIZATION 

AFOSR 


\  8c  ADDRESS  (City,  State,  and  ZIP  Code) 


8b.  OFFICE  SYMBOL 
(If  applicable) 

NL 


Form  Approved 
OMB  No.  0704-0188 


DISTRIBUTION /AVAILABILITY  OF  REPORT 
pproved  for  public  release; 
stribution  unlimited 


MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

A'FOSR-TR'  9  5-0175  _ 


7a.  NAME  OF  MONITORING  ORGANIZATION 

Air  Force  Office  of  Scientific  Research 


7b.  ADDRESS  {City,  State,  and  ZIP  Code) 
Building  410 

Bolling  AFB  D  C  20332-6448 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

F49620-92- J-0510  _ 


10.  SOURCE  OF  FUNDING  NUMBERS 


1 


PROGRAM 
ELEMENT  NO. 


PROJECT 

NO. 


Building  410 

Bolling  AFB  D  C  20332-6448 


11.  TITLE  (Include  Security  Classification) 

Eurasian  Seismic  Surveillance  -  2D  FD  Seismic  Synthetics  and 

Event  Discrimination  _ -  _ 


12.  PERSONAL  AUTHOR(S) 

Bent  0.  Ruud,  Eystein  S.  Husebye  _ 


13a.  TYPE  OF  REPORT  13b.  TIME  COVERED  Il4.  DATE  OF  REPORT  (Year,  Month,  Day)  15.  PAGE  COUNT 

Final  Technical  Report  from  i_  Oct  92  to^Q  SeP  1995  February  23 

16.  SUPPLEMENTARY  NOTATION 


WORK  UNIT 
ACCESSION  NO. 


COSATI  CODES  l  18.  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

SUB-GROUP  2D  FD  synthetics  including  topography,  inversion  for  crustal 

structure,  real-time  magnitude  estimation,  seismic  event 
classification 


19.  ABSTRACT  { Continue  on  reverse  if  necessary  and  identify  by  block  number)  1 

Most  efforts  have  been  devotes  to  2D  FD  synthetics  for  simulating  seismic  wave 
propagation  in  a  complex  lithosphere.  The  FD  scheme  used  is  tied  to  the  numerical 
solution  of  the  el  as  to  dynamic  wave  equation.  We  have  demonstrated  the  usefulness  of 
this  approach  by  computing  synthetic  seismograms  for  a  set  of  progressively  complex 
lithosphere  models.  In  practice,  model  set-up  for  FD  computations  have  proved 
cumbersome  so  we  have  elaborated  on  automating  structural  model  generation.  The 
computed  synthetics  are  SEG-Y  formatted  and  thus  are  easily  processed  using  standard 
software  packages  like  the  ProMAX. 

As  part  of  the  2D  FD  studies,  we  have  critically  examined  optimal  design  of 
boundary  conditions  in  order  to  reduce  edge  reflection  artefacts.  Our  synthetic 
wavefield  scheme  recently  extended  to  3D  modelling  can  also  handel  surface  topography 

>  (continues  on  next  page) 


20.  DISTRIBUTION /AVAILABILITY  OF/A8STRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

□  unclassified/unlimited  121  same  as  rpt.  □  dtic  USERS  Unclassified 


.  22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

Dr.  Stanley  K.  Dickinson 


DO  Form  1473,  JUN  86 


Previous  editions  are  obsolete 


21.  ABSTRACT  SECURITY  CLASSIFICATION 

Unclassified  _ 


22b.  TELEPHONE  (Include  Area  Code)  22 c.  OFFICE  SYMBOL 

(202)  767-4964  AFOSR/NL _ 


bsoiete.  SECURITY  CLASSIFICATION  OF  THIS  PAGE 

.  ^ Unclassified 


19950328  10 


UNCLASSIFIED 


In  the  latter  case,  we  have  been  able  to  simulate  in  3D  Rg-scattering  observations  as 
observed  at  the  NORESS  array.  Crustal  wave  propagation  is  important  in  a  monitoring 
context  so  we  have  studied  in  a  realistic  manner  the  relevance  of  overly  detailed 
crustal  models  frequently  published  in  the  scientific  litteratur.  These  experiments 
were  tied  to  synthetic  wavefields  generated  both  for  homogeneous  and  inhomogeneous 
(3  percent  pertubation) ,  and  then  trying  to  recover  the  original  models  through 
tau-p  related  analysis  and  inversion.  Even  for  oversimplified  homogeneous  crustal 
models  (no  P-  and  S-coda  excitation)  we  were  not  able  to  recover  the  original 
crustal  model  in  detail.  We  concluded  that  a  crust  with  linearized  velocity 
increases  in  combination  with  2-4  percent  velocity  pertubations  is  a  better 
approximation  to  real  Earth  structures  than  the  often  published,  homogeneous 
multilayered  models. 

Event  magnitude  is  an  versatile  parameter  in  many  seismo logical  studies  not  at 
least  in  the  context  of  nuclear  test  ban  monitoring.  For  simulating  network  signal 
detection  capabilities  the  noise  level  at  individual  stations  are  needed  in  near 
real  time  which  in  practice  appears  to  be  difficult  to  obtain.  Ideally,  we  would 
prefer  to  have  such  noise  level  estimates  available  on-line  to  measure  the  * on-line 
threshold  monitoring  performance  and  likewise  event  magnitudes .  To  achieve  this  we 
have  been  exploring  the  usefulness  of  a  random  vibration  theory 
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1  Objectives 

•  Foster  a  better  understanding  of  seismic  wave  propagation  in  complex  lithospheric  media. 

•  Principal  modeling  tool  is  2D  finite  difference  (FD)  synthetics  for  models  also  including 
surface  topography. 

•  Corroborating  synthetics  through  analysis  of  array  and  network  data. 

•  Use  ’synthetic’  knowledge  in  design  of  event  classification  schemes. 

•  Seismic  network  performance. 

Note,  a  basic  element  in  research  is  curiosity  which  in  practice  will  imply  a  certain  flexibility 
regarding  adherence  to  objectives. 


2  Report  summary 

Most  efforts  have  been  devotes  to  2D  FD  synthetics  for  simulating  seismic  wave  propagation 
in  a  complex  lithosphere.  The  FD  scheme  used  is  tied  to  the  numerical  solution  of  the  elastody- 
namic  wave  equation.  Spatial  partial  differentiation  is  achieved  through  cost-optimized,  dispersion- 
bounded,  high-order  finite  difference  operators  on  a  staggered  grid.  For  time  stepping,  a  leap-frog 
technique  is  used.  For  the  numerical  dispersion  relations,  the  stability  limit  and  bandwidth  intro¬ 
duced  by  the  discretization,  the  reader  is  referred  to  Sguazzero,  Kindeland  and  Kamel  (1989). 

Recently  we  have  demonstrated  the  usefulness  of  this  approach  by  computing  synthetic  seis¬ 
mograms  for  a  set  of  progressively  complex  lithosphere  models  (Hestholm,  Husebye  and  Ruud, 
1994).  In  practice,  model  set-up  for  FD  computations  have  proved  cumbersome  so  we  have  elab¬ 
orated  on  automating  structural  model  generation.  The  new  approach  now  working  satisfactorily 
is  tied  to  using  a  scanner  for  reading  say  a  published  multilayered  crustal  model.  For  any  model 
point,  proper  layer  assignment  is  performed  automatically  and  hence  a  corresponding  specifica¬ 
tion  of  phvsical  parameters  like  velocity  and  density  is  obtained.  The  computed  synthetics  are 
SEG-Y  formatted  and  thus  are  easily  processed  using  standard  software  packages  like  the  popular 
ProMAX. 

The  method  of  r-sum  inversion  provide  a  convenient  way  to  estimate  velocity-depth  models 
from  seismic  profiling  data  via  r-p  transformation.  Two  basic  assumptions  must  be  valid:  (i)  the 
velocity  must  increase  monotonically  with  depth,  and  (ii)  the  medium  must  be  laterally  homoge¬ 
neous.  In  a  new  paper  we  test  the  sensitivity  of  the  inversion  method  for  small  deviations  from 
the  second  of  these  assumptions  by  inverting  synthetic  seismograms  computed  for  two  models,  one 
laterally  homogeneous  and  the  other  perturbed  by  a  2-dimensional  random  field.  The  synthetics 
were  generated  by  finite  difference  modeling  of  the  elasto-dynamic  wave  equation.  Also  the  gen¬ 
eral  resolution  of  the  method  for  the  actual  frequency  bandwidth  (1-5  Hz)  and  the  experimental 
layout  (1  shot,  sensors  at  0.4  km  intervals  out  to  320  km  offset)  could  be  attained.  The  results 
showed  that  for  the  perturbations  used  (von  Karman  correlation  function  of  order  0.3,  correlation 
distances  of  10  km  and  2.5  km  in  the  horizontal  and  vertical  directions  respectively,  and  an  RMS 
of  3%  )  the  model  could  be  recovered  with  about  the  same  accuracy  whether  perturbed  or  not. 
The  general  resolution  in  velocity  was  about  0.2  km/s  all  the  way  down  to  the  Moho.  The  results 
were  somewhat  surprising  in  view  of  the  rather  strong  scattering  produced  by  the  perturbations, 
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wiping  out  most  of  the  phases  observable  for  the  unperturbed  model.  The  success  of  the  inversion 
method  in  recovering  the  large  scale  structure  of  the  velocity-depth  function  is  attributed  to  the 
fact  that  it  depends  only  on  the  strongest  post-critically  reflected  phases.  However,  fine  scale 
layering  of  the  models  was  unrecoverable  even  for  the  unperturbed  model. 

An  important  motivation  for  this  particular  study  was  is  a  disbelief  in  the  realism  of  the 
many  finely  layered  crustal  models  derived  from  large-  scale  wide-angle  reflection  and  refraction 
profiling  experiments.  Simply,  the  claimed  crustal  models  hardly  have  a  counterpart  vis-a-vis  the 
field  observations  at  hand.  Our  disbelief  in  such  models  have  been  further  strengthened  during  the 
experiments  described  here  -  under  ideal  simulated  observational  conditions  like  no  background 
noise,  dense  spatial  wavefield  sampling  and  using  typical  crustal  layer  velocities  -  still  we  are  unable 
to  retrieve  the  finer  details  of  the  a  priori  known  crustal  model.  Adding  velocity  perturbations,  a 
prerequisite  for  generating  realistic  synthetics  in  terms  of  P-  and  S-coda  excitation,  crustal  layering 
details  are  definitely  non-retrievable.  To  summerize;  simple  crustal  models  in  terms  of  linearized 
velocity  increases  with  depth  and  in  combination  with  structural  perturbations  of  2-  4  per  cent 
appear  to  match  real  Earth  structures  far  better  than  elaborated,  multilayered  crustal  models. 

Event  magnitude  is  an  versatile  parameter  in  many  seismological  studies  not  at  least  in  the 
context  of  nuclear  test  ban  monitoring.  For  simulating  network  signal  detection  capabilities  the 
noise  level  at  individual  stations  are  needed  in  near  real  time  which  in  practice  appears  to  be 
difficult  to  obtain.  Ideally,  we  would  prefer  to  have  such  noise  level  estimates  available  on-line 
to  measure  the  ’on-line  threshold  monitoring5  performance  (Ringdal  and  Kvaerna,  1989).  To 
achieve  this  we  have  been  exploring  the  usefulness  of  a  random  vibration  theory  formula,  namely 
(Cartwright  and  Longnet-Higgins,  1956). 

Amax  =  f(N)A  rms  (i) 

where  Amax  is  max  amplitude;  f(N)  —  iV(2  In  AQ1/2  where  N  is  the  number  of  extrema  with 
the  RMS- window  and  Arms  is  amplitude  root-mean-square.  Most  signal  detections  in  use  are 
of  the  STA/LTA  type  where  the  short-term-average  (STA)  is  in  fact  a  root-mean-square  (RMS) 
trace  estimate.  Initial  test  results  from  real  recordings  indicate  that  the  above  formula  is  valid; 
error  estimates  amount  to  ±0.1  magnitude  unit.  In  other  word,  on-line  magnitude  estimation 
should  be  feasible.  Extensive  analysis  of  explosion  and  earthquake  recordings  by  the  Norwegian 
Seismological  Network  indicated  that  this  is  valid  for  practical  seismological  purposes.  On  this 
basis,  we  undertook  local  Ml  (Lg)  magnitude  estimations  for  30  events;  satisfactory  results  were 
obtained  as  our  novel  Ml  determinations  were  very  similar  to  those  reported  by  NORSAR  (Mendi 
and  Husebye,  1994). 

Local  magnitude  estimation  is  an  important  issue  in  the  seismic  monitoring  context,  but  in 
practise  often  somewhat  problematic  due  to  lack  of  proper  empirical  signal  frequency-distance 
correction  curves.  To  overcome  this  kind  of  problems,  we  constructed  ‘theoretical’  frequency- 
distance  curves  for  Scandinavia  using  attenuation  and  geometrical  spreading  constants  as  pub¬ 
lished  by  Sereno  et  al  (1988).  Good  agreement  was  obtained  with  conventional  correction  curves 
thus  demonstrating  that  local  magnitude  estimation  procedures  are  transportable  conditioned  on 
available  geometric  spreading  and  attenuation  constants. 

Finally,  a  comprehensive  study  aimed  at  seismic  event  classification  at  local  and  regional  dis¬ 
tances  was  considered.  Preparatory  efforts  here  comprise  software  development  and  selection  of 
a  comprehensive  data  base;  we  opted  for  the  NORESS  and  ARCESS  events  used  by  Lacoss  et 
al.  (1992).  Regrettably,  we  cannot  report  much  progress  on  this  issue  under  this  contract.  What 
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we  have  done  so  far  was  to  compute  the  P/Lg  amplitude  ratios  for  the  30  events  analysed  in 
the  Mendi  and  Husebye  (1994)  study,  also  in  an  automated  manner.  The  results  obtained  were 
not  very  positive  as  regards  this  particular  discriminant  as  there  were  no  clear  distinction  be¬ 
tween  the  earthquake  and  explosion  populations.  It  is  probably  several  reasons  for  this;  firstly, 
source-receiver  paths  are  often  offshore  in  a  continental  crustal  regime  with  many  thick  sedimen¬ 
tary  basins;  secondly,  the  'ground  truth’  aspect  of  the  Bergen  University  bulletin  is  poor,  as  the 
earthquakes  and  explosions  are  frequently  intermixed.  Close  inspection  of  our  results  indicate  that 
the  very  calculations  of  P/Lg  discriminant  parameters  via  the  Arms  relation  given  in  eq.(l)  would 
provide  in  near  real  time  information  bearing  on  the  source  identification. 


2.1  Outline  of  the  report 

In  the  next  section  of  this  Final  Technical  Report  under  Contract  AFOSR  Grant  F49620-92-5- 
0510  we  describe  the  most  recent  research  accomplished  concerning  modeling  of  scattering  from 
topographic  relief  near  the  NORESS  array.  This  work  is  performed  with  a  new  3D  elasto-dynamic 
finite  difference  code  allowing  for  topography  at  the  free  surface.  The  theoretical  details  of  this 
novel  method  is  given  in  Appendix  A.  We  plan  to  continue  to  work  on  this  interesting  modeling 
experiments  under  the  AFOSR  Grant  F49620-94-1-0278  at  the  University  of  Bergen.  In  Appendix 
B  we  include  a  new  paper  by  Simone  and  Hestholm  (1995)  concerning  the  absorbing  boundary 
method  used  in  our  finite  difference  programs.  The  complete  text  of  the  Ruud  et  al  (1995)  paper 
of  tau-p  inversion  of  synthetic  2D  FD  refraction  and  wide-angle  reflection  profiles  is  given  in 
Appendix  C,  and  in  Appendix  D  we  include  the  already  published  paper  by  Mendi  and  Husebye 
(1994)  on  near  real  time  magnitude  estimation. 


References 

Cartwright,  D.E.,  and  Longuet-Higgins,  M.S.,  1956.  The  statistical  distribution  of  the  maxima  of 
a  random  function,  Proc.  Roy.  Soc.  London,  Ser.  A237 ,  212-223. 

Hestholm,  S.O.,  Husebye,  E.S.,  and  Ruud,  B.O,  1994.  Seismic  wave  propagation  in  complex  crust 
-  upper  mantle  media  using  2D  finite  difference  synthetics,  Geophys.  J.  Ini.,  118,  643-670. 

Hestholm,  S.O.  and  Ruud,  B.O,  1994.  2D  finite  difference  elastic  wave  modeling  including  surface 
topography,  Geophyscal  Prospecting,  42,  371-390. 

Lacoss,  R.,  Cunningham,  R.,  Curtis,  S.,  and  Seibert,  M.,  1991.  Artificial  neural  networks  for 
seismic  interpretation,  Technical  Report,  ESD-TR-91-170,  MIT  Lincoln  Laboratory,  Lexington, 
MA. 

Mendi,  C.D.  and  Husebye,  E.S,  1994.  Near  real  time  seismic  magnitude  and  moment  measure¬ 
ments,  Annali  di  Geofisica ,  37,  365-382. 

Ringdal,  F.,  and  Kvaerna,  T.,  1989.  A  multi-channel  processing  approach  to  real  time  network 
detection,  and  threshold  monitoring,  Bull.  Seism.  Soc.  Am.,  79,  1927-1940. 

Ruud,  B.O,  Grundt,  Y.,  and  Husebye,  E.S.,  1995.  A  comparison  of  tau-sum  velocity-depth  inver¬ 
sion  for  a  perturbed  and  unperturbed  model,  m/s  prepared  for  publication. 


5 


Sereno,  T.J.,  Bratt,  S.R.,  and  Bache,  T.C.,  1988.  Simultaneous  inversion  of  regional  wave  spectra 
for  attenuation  and  seismic  moment  in  Scandinavia,  J.  Geophys.  Res.,  93,  2019-2035. 

Sguazzero,  P.,  Kindeland,  M.,  and  Kamel,  A.,  1989.  Dispersion-Bounded  Numerical  Integration  of 
the  Elastodynamic  Equations,  Proceedings  of  the  ICOSAHOM  Conference ,  Como,  Italy,  June, 
1989 ,  pp.  165-172.  North  Holland  Publishing  Company,  The  Netherlands. 

Simone,  A.  and  Hestholm,  S.,  1995.  Optimal  absorbing  boundary  condition  in  2D  finite  differ¬ 
ence  elastic  wave  modeling.  Empirical  study  and  comparison  with  exponential  damping,  m/s 
prepared  for  publication. 


3  3D  FD  modeling  of  scattering  from  topographic  relief 

As  in  the  2D  case,  we  expect  that  inclusion  of  topography  at  the  free  surface  will  lead  to  improved 
modeling  of  upper  crust  propagation  effects.  Of  particular  interest  here  is  the  scattering  of  seismic 
energy  from  body  waves  (P  or  S)  into  surface  waves  (mainly  Rg).  This  phenomena  have  been 
extensively  studied  and  documented  using  data  from  the  NORESS  array  in  south-eastern  Norway 
(Bannister  et  al,  1990;  Gupta  et  al,  1993;  Hedlin  et  al,  1991,  1994).  For  our  3D  FD  modeling 
experiments  we  have  obtained  digital  elevation  data  for  an  area  of  100  x  100  km  centered  on  the 
NORESS  array.  Due  to  the  huge  computer  memory  requirements  of  3D  FD  methods,  we  have  so 
far  been  restricted  to  a  model  of  size  40  x  40  x  35  km  with  0.2  km  sampling.  The  topography  in 
the  model  area  is  shown  in  Fig.  1. 


Figure  1:  Topography  in  an  40  x  40  km  area  centered  at  the  NORESS  array.  Horizontal  axes 
are  labeled  in  kilometers.  Elevations  are  in  meters  above  mean  sea  level.  The  dark  area  in  the 
south-west  corner  of  the  figure  is  part  of  Lake  Mjpsa  (123  m  above  sea  level). 
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The  theoretical  background  for  the  implementation  of  3D  FD  code  with  free  surface  topography 
is  given  in  Appendix  A.  This  method  is  a  plane  extension  from  2D  to  3D  of  the  method  described 
by  Hestholm  and  Ruud  (1994).  The  absorbing  boundary  condition  used  along  the  sides  and  bottom 
of  the  model  is  the  socailed  ‘exponential  damping’  technique.  A  description  of  this  method,  and 
comparison  with  another  method,  is  given  in  a  paper  by  Simone  and  Hestholm  (submitted  for 
publication)  included  in  Appendix  B. 

In  all  the  examples  shown  here  the  incoming  wave  is  a  vertical  incident  plane  P-wave  simulating 
a  teleseismic  short  period  P-phase.  The  center  frequency  of  the  Ricker  wavelet  is  2.5  Hz,  the  P~ 
wave  velocity  of  the  homogeneous  medium  is  6.0  km/s  and  the  Poisson’s  ratio  is  0.25.  The  main 
problem  encountered  in  our  test  runs  is  instabilities  which  start  near  the  surface  in  an  area  of 
steep  and  rough  topography.  The  instabilities  appear  as  surface  waves  of  exponentially  increasing 
energy  which  slowly  propagate  out  from  the  starting  point.  In  order  to  quantify  and  locate  areas 
which  cause  instabilities  we  have  computed  different  functions  depending  on  first  and  second  order 
derivatives  of  the  topography  function  2  =  f(x ,  y).  First,  we  computed  the  slope,  i.e.,  the  length 
of  the  topography  gradient  vector  as 


|V/(*,y)| 


Although  the  gradient  of  the  topography  is  useful  to  define  areas  of  strong  scattering,  the  insta¬ 
bilities  of  the  FD  scheme  is  more  likely  related  to  its  roughness.  As  a  measure  here  we  computed 
the  Frobenius  norm  of  the  Hessian  matrix  of  f(x,  y),  i.e., 


\\H(x,y)\\F 


In  Fig.  2  we  show  both  \Vf(x,y)\  and  ||iJ(x,  2/)||f  for  the  topography  in  Fig.  1.  The  derivatives 
were  computed  via  Fourier  transforms.  In  Fig.  2  there  is  a  ‘wall’  of  high  gradient  about  10  km 
east  and  2  km  north  of  the  NORESS  array  (‘Bronkeberget’,  Bannister  et  al,  1990).  This  area  also 
have  the  roughest  topography.  The  area  close  to  Lake  Mjpsa  in  the  south-west  corner  is  within  the 
absorbing  boundaries  (thickness  4  km)  of  our  FD  modeling  scheme  and  was  therefore  not  properly 
modeled  here. 


3.1  Results  and  discussion 

The  first  modeling  attempt  with  unsmoothed  real  topography  gave  very  unstable  results  (infinite 
amplitude)  within  a  second  after  the  wave  front  reached  the  surface.  In  our  next  attempt  the 
elevations  were  first  multiplied  by  0.5  in  order  to  reduce  the  topography.  The  results  were  stable 
and  this  is  also  consistent  with  2D  FD  modeling  results.  In  the  latter  case  an  east-west  profile 
running  through  NORESS  was  unstable  with  real  topography  but  stable  with  half  elevations. 
Snapshots  of  the  wavefield  (vertical  component  of  the  particle  velocity)  at  the  free  surface  is 
shown  in  Fig.  3. 

A  dominant  feature  of  these  snapshots  is  the  low-frequency /long- wavelength  artificial  reflections 
from  the  absorbing  boundaries.  Although  we  consider  the  exponential  damping  technique  to  be 
the  most  efficient  absorbing  boundary  method  we  have  tested,  the  case  of  an  wave  incident  at 
an  angle  close  to  90°  with  the  boundary  is  the  most  difficult  one.  In  our  snapshots  the  artificial 
reflections  are  seen  to  generate  surface  waves  starting  at  the  absorbing  boundaries,  4  km  from  the 
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Figure  3:  Snapshots  of  the  vertical  component  of  the  particle  velocity  at  the  free  surface.  The 
time  in  the  upper  left  corner  of  each  snapshot  is  the  time  from  when  the  plane  vertical  incident 
P-wave  reached  the  surface.  The  first  snapshot  (at  0.0  sec)  looks  different  from  the  others  because 
the  average  amplitude  differs  from  zero.  The  straight  wavefronts  parallel  to  the  model  boundaries 
are  surface  waves  generated  artificially  from  the  absorbing  boundaries.  Note  the  strong  scattering 
source  located  at  Bronkeberget  (2  km  N,  10  km  E)  which  is  also  highly  prominent  in  real  record 
analysis  (Bannister  et  al,  1990). 
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edges  of  the  model.  Fortunately,  the  artificial  reflections  are  frequency  dependent,  with  longer 
wavelengths  than  most  of  the  surface  waves  scattered  by  the  topography,  and  can  therefore  be 
removed  by  image  processing  techniques.  This  have  been  done  by  applying  a  high-pass  filter  in 
2D  wavenumber  domain  which  removes  waves  with  wavelengths  larger  than  2.4  km  and  leaves 
waves  with  wavelengths  less  than  1.2  km  —  a  cosine  taper  is  used  in  the  transition  band.  The 
snapshots  after  wavenumber  filtering  is  shown  in  Fig.  4.  As  seen  from  this  figure  the  scattered 
surface  waves  appear  to  radiate  out  from  secondary  point  sources  which  coincide  with  areas  of 
high  topographic  gradients  (Fig.  1).  The  dominant  scattering  points  are  along  the  steep  valley 
side  east  of  ‘Bronkeberget’  about  10  km  east  of  NORESS. 

As  the  instabilities  in  our  3D  FD  modeling  is  likely  due  to  the  roughness  of  the  free  surface,  we 
smoothed  the  topography  to  see  if  stable  results  could  be  obtained  without  reducing  the  elevation 
of  the  large  scale  topographic  features.  A  low-pass  version  of  the  2D  wavenumber  filter  was  used 
to  smooth  the  topography.  Several  models  were  tested  and  it  appeared  that  a  rather  high  degree 
of  smoothing  was  necessary  to  obtain  stable  results.  The  roughest  topography  model  that  seems 
to  be  stable  was  generated  by  removing  all  wavelengths  shorter  than  1.2  km  —  the  cosine  tapered 
transition  band  started  at  2.4  km.  The  smoothed  topography  model  is  shown  in  Fig.  5  and  its 
derivatives  in  Fig.  6.  The  snapshots  computed  for  the  new  smoothed  model  are  given  in  Fig.  7 
and  the  filtered  versions  in  Fig.  8.  In  this  case  the  scattering  points  are  also  along  a  line  about 
10  km  east  of  NORESS,  but  now  extending  further  south  compared  to  the  unsmoothed  (but 
reduced)  topography.  This  result  might  be  expected  as  also  the  gradient  map  (Fig.  6)  exhibits 
this  southward  extension.  However,  more  surprising  is  the  fact  that  the  scattered  waves  propagate 
mainly  in  the  eastward  direction  -  the  unsmoothed  model  produced  no  distinct  radiation  pattern. 

The  amplitudes  of  the  scattered  waves  were  about  the  same  in  the  two  modeling  experiments 
shown  here.  The  gradients  were  also  of  the  same  size  when  we  take  into  account  that  the  values 
in  Fig.  2  must  be  multiplied  by  0.5  to  be  consistent  with  the  model  actually  used.  However, 
the  Frobenius  norm  of  the  Hessian  matrix  (computed  from  the  second  order  derivatives  of  the 
topography)  is  much  smaller  for  the  smoothed  model.  Also  for  other  smoothed  models  (not  as 
smooth  as  the  one  discussed  here)  for  which  the  Frobenius  norm  was  much  smaller  than  for  the 
reduced  but  unsmoothed  model,  the  modeling  gave  instable  results.  It  thus  appear  that  the 
roughness  of  the  topography  is  not  the  only  reason  for  the  instability  or  that  another  measure 
of  the  roughness  is  required.  More  work,  including  modeling  with  simpler  ‘synthetic’  types  of 
topography,  is  necessary  in  order  to  understand  the  origin  of  the  instability  and  possibly  eliminate 
the  problem. 

3.2  Concluding  remarks 

Even  for  small  scale  experiments  3D  FD  synthetic  wavefield  analysis  provide  as  demonstrated 
here  improved  insight  and  a  better  understanding  of  surface  scattering  phenomenon.  Particular 
gladdening  is  that  the  strong  Rg  scattering  hill  Bronkeberget,  identified  as  such  by  Bannister  et  al. 
(1990)  through  analysis  of  NORESS  recordings,  can  be  realistically  synthesized.  Another  inter¬ 
esting  feature  is  that  for  certain  spatial  low-frequency  variants  of  the  local  NORESS  topography 
the  wave  propagation  characteristics  could  change  abruptly  -  the  eastward  only  propagation  of 
Rg-waves  excited  by  secondary  sources  in  the  hilly  Bronkeberget  area.  Such  wavefield  features 
are  reminiscent  of  NORESS  and  GERESS  recordings;  at  some  sensors  the  Rg-phase  is  prominent 
while  hardly  visible  at  nearby  sensors  less  than  a  kilometer  away.  In  other  w'ords,  also  small  scale 
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Figure  4:  Wavenumber  filtered  versions  of  the  snapshots  in  Fig.  3.  The  surface  waves  artificially 
generated  at  the  absorbing  boundaries  are  removed  by  high-pass  filtering  above  2.4  km  wave¬ 
lengths.  The  remaining  surface  waves  are  generated  by  scattering  of  the  topographic  relief  (see 
Fig.  1).  As  in  Fig.  3,  the  Bronkeberget  scattering  hill  remains  dominant. 
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Figure  5:  Smoothed  topography  in  an  40  x  40  km  area  centered  at  the  NORESS  array  obtained 
by  2D  wavenumber  filtering  -  a  low-pass  filter  removing  all  wavelengths  less  than  1.2  km  was  used. 
Horizontal  lengths  are  in  kilometers  while  elevations  are  in  meters. 

crustal  features  may  contribute  to  crustal  blocking  effects  as  often  observed  for  Rg-  and  Lg-phases. 
Many  of  such  wave  field  phenomena  are  out  of  range  for  2D  FD  wave  synthetic  experiments  hence 
our  emphasize  on  continuous  AFOSR  contracted  research  efforts  on  3D  synthetic  simulation  of 
crustal  wave  field  propagation. 


Bent  0.  Ruud 


Stig  Hestholm 


Eystein  S .  Husebye 
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Figure  7:  Snapshots  of  the  vertical  component  of  the  particle  velocity  at  the  free  surface  computed 
for  the  smoothed  topography  model  (Fig.  5).  The  first  snapshot  (at  0.0  sec)  looks  different  from 
the  others  because  the  average  amplitude  differs  from  zero.  The  straight  wavefronts  parallel  to  the 
model  boundaries  are  surface  waves  generated  artificially  from  the  absorbing  boundaries.  Note 
that  the  Bronkeberget  scattering  area  (2  km  N,  10  km  E)  produce  strongly  directional  dependent 
(eastward  propagating)  Rg  waves. 
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Figure  8:  Wavenumber  filtered  versions  of  the  snapshots  in  Fig.  7.  The  surface  waves  artifi¬ 
cially  generated  at  the  absorbing  boundaries  are  removed  by  high-pass  filtering  above  2.4  km 
wavelengths.  The  remaining  surface  waves  are  mainly  generated  by  scattering  of  the  smoothed 
topographic  relief  (Fig.  5). 
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Appendix  A 


Theoretical  Background  for 
3D  Finite  Difference  Elastic  Wave  Modeling 
including  Surface  Topography 


Introduction 

In  this  work  explicit  3D  boundary  conditions  for  a  free  surface  topography  are  deduced. 
The  basis  is  the  vanishing  stress  condition  for  a  free  surface.  A  3D  curved  grid  which  is 
stretched  in  the  vertical  direction  is  used  to  adapt  to  the  surface  topography,  i.  e.  the  top 
surface  of  the  grid  coincides  with  the  surface  topography.  A  coordinate  transform  is  used 
to  transfer  the  elastic,  isotropic  wave  equations  from  the  curved  to  a  rectangular  grid  in 
which  the  numerical  computations  are  done. 

At  the  topography  surface,  the  velocity  boundary  conditions  for  a  free  surface  are 
implemented  into  a  local,  rotated  system  at  each  point  on  the  surface.  Each  of  these 
systems  has  its  vertical  coordinate  direction  coinciding  with  the  normal  vector  of  the 
surface  at  the  given  point.  The  velocity  boundary  conditions  are  subsequently  rotated 
back  to  the  rectangular  system.  Once  the  boundary  conditions  are  given  in  this  system, 
the  numerical  computations  can  be  performed. 

In  the  following  paragraphs,  the  3D  equivalents  to  the  2D  equations  given  in  Hestholm 
and  Ruud  (Geophysical  Prospecting,  1994,  42:371-390)  will  be  stated.  Brief  descriptions 
of  the  deductions  will  also  be  given.  For  each  set  of  equations,  the  corresponding  equation 
numbers  for  the  2D  case  in  Hestholm  and  Ruud  (1994)  will  be  stated.  In  each  instance,  it 
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is  possible  to  verify  the  coincidence  between  the  3D  and  the  2D  case  by  eliminating  either 
of  the  horizontal  directions  in  3D.  Also  stated  is  the  coincidence  between  the  conditions 
for  the  surface  topography  applied  to  a  plane  surface,  and  the  plane  surface  conditions  in 
3D. 

Elastic  Wave  Modeling  Formulation 

The  linear  transformation  from  the  rectangular  (£,  r,  ^-system  to  the  stretched  (x,y,z)~ 
system,  where  both  systems  have  positive  direction  upwards  along  the  vertical  coordinate, 
can  be  written  in  3D  as  (Hestholm  and  Ruud,  1994,  eqs.  (4)  and  (5)) 

«(£.r,»7)  =  £,  (!') 

y($>T,rj)  =  T,  (2') 

*(£»  D  v)  =  r),  (3') 

T}max 


where  zo(£)r)  is  the  topography  function. 

The  expressions  for  the  partial  derivatives  corresponding  to  eqs.  (8) — (10)  in  Hestholm 
and  Ruud  (1994),  which  are  needed  in  the  medium  equations,  are  found  to  be 
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We  expand  eqs.  (l)-(3)  in  Hestholm  and  Ruud  (1994)  by  the  chain  rule  and  put  the 
dimension  J  =  3.  Then  we  get  the  equations  for  the  wave  propagation  inside  the  medium 
corresponding  to  eqs.  (11)— (15)  in  the  Hestholm  and  Ruud  (1994), 

du  d<7 xx  ,  t  ,  d<JXy  n/V  ^  d(Txy  rt(e  ,  t  /fti\ 

P^7  =  —  +  A^,r,r])-^-  +  -^  +  B^T,r1)—+C(tT)—  +  fx,  (9  ) 


at  at 


fTt  =  it + A((’ r- ^ + ^ + m T' + C(?’  T)^r + fr •  im 

=  w + Ai(’- r’ n)S-w + ^ + B(f  ■  ■ r' 1 ^ w + • ■ r)d~w + ■ (1 1 r) 

Tf  =  (A  +  2,)  (|  +  A({,r,,)^)+A(|^  +  B(?,  r,  ,)|  +  C«,  r)£)  ,  (12') 

=  A(w  +  '4K,r',')^  +  c'K’r)^)+(A  +  2,‘)(^  +  SK’r’’')^)’ (13<) 
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^  =  ,(^  +  ^.r„)^  +  C«,^),  (16') 

^=,(^  +  B(C^)|+C«,r)|).  (17') 


where  p  is  the  density  and  A  and  p  are  Lame’s  parameters.  fX)  fy  and  fz  are  the  components 
of  the  body  forces,  and  u,  v  and  w  are  the  particle  velocity  components.  crXX)  <7yy,  azz. 


*  xy 3  u  xz 


crxz  and  avz  are  the  stress  components. 


Free  Surface  Boundary  Conditions 

The  3D  free  boundary  conditions  on  the  velocities  at  a  locally  horizontal  surface  (or  where 
the  z— axis  coincides  with  the  local  normal  vector  of  the  surface)  resulting  from  the  vanishing 
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stress  condition  can  be  written 


A  fdu  dv 

A  +  2/i  \dx  dy 


with  x  and  y  the  horizontal  coordinates  and  z  the  vertical  coordinate.  These  equations 
correspond  to  eqs.  (16)  and  (17)  in  Hestholm  and  Ruud  (1994). 

We  want  to  apply  these  conditions  to  a  topography  surface.  At  each  surface  point , 
they  then  have  to  be  imposed  into  a  local  coordinate  system  (xf ,  y',  z')  in  which  the  z'-axis 
coincides  with  the  local  normal  vector  of  the  surface.  In  this  local  system  we  impose  the 
conditions  (18/)-(20/).  Once  this  is  done,  they  have  to  be  rotated  back  to  the  rectangular 
system  (£,r,  rf)  in  order  to  do  the  numerical  computations.  This  rotation  is  expressed  by 

v  =  Av  ',  (21') 


where  v  and  v  '  are  the  particle  velocity  vectors  in  the  (£,  r,  77)-  and  the  (xl ,  yf ,  z^-systems 
respectively.  A  is  the  rotation  matrix.  This  was  found  by  geometrical  considerations,  and 
is  given  by 


cos  9  0  sin  9 

-  sin  9  sin  <j>  cos  6  cos  9  sin  0 
■  sin  9  cos  (f)  —  sin  0  cos  9  cos  <p 


9  is  the  rotation  angle  between  the  £-axis  and  the  local  z'-axis  in  the  (£,  ?7)-plane.  6  is 
the  rotation  angle  between  the  r-axis  and  the  local  y'-axis  in  the  local  (?/,  z/)-plane.  We 
also  have  the  relations  tan0  =  dzo(^)r)/d^  and  tan  6  —  dzo(£,  r)fdr. 

Now  the  calculations  of  the  rotation  from  the  (x* ,  y\  z^-system  back  to  the  (£,  r,  r] )- 
system  are  performed  in  a  completely  analogue  way  as  is  done  in  Appendix  3  of  Hestholm 


and  Ruud  (1994)  for  the  2D  case.  We  arrive  at  the  3D  boundary  conditions  for  a  free 
surface  topography  given  in  the  computational  (£,  r,  7?)-system  by 


(29') 

(30') 

(31') 


Numerical  Discretization 

For  the  numerical  discretization,  the  reader  is  referred  to  the  corresponding  paragraph 
in  Hestholm  and  Ruud  (1994).  The  same  spatial  and  time  discretization  procedures  are 
used  as  in  the  2D  case.  In  order  to  get  explicit  expressions  for  each  physical  variable  at 
each  time  step,  we  still  have  to  stagger  the  vertical  velocity  component  w  one  half  grid 
length  downwards.  Corresponding  numerical  definitions  have  to  be  implemented.  The 
3D  boundary  conditions  eqs.  (23')-(25')  are  discretized  by  second  order,  central,  staggered 
finite  difference  operators.  Below  the  free  surface,  the  operator  order  is  gradually  increased 
with  depth  up  to  8th  order.  Generally,  u  is  staggered  one  half  grid  length  in  the  positive 
^-direction,  v  is  staggered  one  half  grid  length  in  the  positive  r-direction,  while  w  is 
staggered  one  half  grid  length  in  the  negative  ^-direction  (downwards). 

To  find  the  velocity  components  at  the  surface  topography  from  the  conditions  (23')- 
(25'),  we  solve  directly  a  simultaneous  3X3  system  by  using  Cramer’s  rule.  The  vertical 
grid  level  for  which  the  conditions  are  solved  are  not  consistent  for  all  the  terms  and 
equations,  with  up  to  a  grid  length  difference  vertically.  As  for  the  2D  case,  this  seems 
just  to  smooth  the  boundary  to  some  extent,  and  even  be  to  some  advantage  numerically. 
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Appendix  B 


Optimal  Absorbing  Boundary  Condition  in 
2D  Finite  Difference  Elastic  Wave 
Modelling.  Empirical  Study  and 
Comparison  with  Exponential  Damping 

A.  Simone1  and  S.  Hestholm2’3 

1  Istituto  Universitario  Navale,  Via  Amm.  Acton  38,  Naples,  Italy 

2  IBM  Bergen,  0.  Kr&kenes  17,  5042  Fjpsanger,  Norway 
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Introduction 

The  problem  of  artificial  reflections  from  grid  boundaries  in  the  numerical  discretization 
of  elastic  and  acoustic  wave  equations  has  intrigued  geophysicists  for  a  long  time.  Even  if 
modern  computers  have  made  it  possible  to  extend  the  synthetics  over  more  wavelengths, 
equivalent  to  larger  propagation  distances,  efficient  absorption  methods  are  still  needed  in 
order  to  minimize  interference  from  artificial  grid  boundary  reflections.  In  this  study  we 
examine  the  relative  merits  of  a  recently  published  absorption  method,  namely  the  Opti¬ 
mal  Absorbing  Boundary  Condition  (OABC)  designed  by  Peng  and  Tokspz  (1994)  for  2D 
and  3D  acoustic  and  elastic  wave  modelling.  As  a  basis  for  comparison  we  use  the  Expo¬ 
nential  Damping  (ED),  in  which  velocities  and  stresses  are  multiplied  by  progressively 
decreasing  terms  towards  the  edges.  The  experiment  design  is  to  implement  both  the 
OABC  and  the  ED  in  our  2D  finite  difference  elastic  wave  algorithm,  see  Hestholm  and 
Ruud  (1994),  and  then  check  their  respective  performances.  A  disadvantage  with  ED  is 
that  “absorption  strips”  are  involved,  which  means  in  case  of  3D  modelling,  that  a  signifi¬ 
cant  amount  of  additional  computer  storage  is  required.  On  the  other  hand,  initially  we  had 
some  problems  with  the  implementation  of  the  OABC  into  our  2D  finite  difference 
scheme,  because  instabilities  occurred  in  the  synthetic  wavefield  calculations.  Ways  to 
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delay  the  occurrence  of  these  instabilities  will  be  adressed,  in  addition  to  the  relative  mer¬ 
its  of  the  OABC  and  ED  methods. 

2-D  finite  difference  implementation 

The  2D  finite  difference  scheme  employed  solves  the  basic  equations  governing  wave 
propagation  in  an  elastic  isotropic  medium,  using  the  velocity-stress  formulation, 
(Hestholm  and  Ruud,  1994).  The  discretization  of  the  elastodynamic  equations,  with  two 
staggered  numerical  space  differentiators,  has  been  applied  as  in  Levander  (1988).  Details 
of  the  numerical  discretization  are  given  by  Sguazzero  etaC.  (1989).  .Spatial  partial  differ¬ 
entiation  is  achieved  by  finite  difference  operators  on  a  staggered  grid.  The  scheme  used  in 
the  interior  of  the  computational  domain  is  accurate  to  the  8th  order  in  space.  For  time 
stepping  a  leap  frog  technique  is  used. 

The  method  of  OABC  extrapolates  values  on  the  numerical  boundaries  of  a  finite  differ¬ 
ence  domain.  It  is  expressed  as  a  linear  combination  of  wave  fields  at  previous  time  steps 
and/or  interior  grids,  and  it  makes  use  of  the  zeros  and  poles  of  the  reflection  coefficients 
in  the  complex  plane.  Peng  and  Tokspz  (1994)  pointed  out  that  the  purpose  of  the  OABC 
was  to  design  an  absorbing  boundary  condition  that  was  optimal  in  the  sense  that  “maxi¬ 
mum  absorbing  effects  can  be  achieved  with  a  minimum  amount  of  computation  and  stor¬ 
age,  and  is  easily  implemented”.  Approaching  the  boundaries,  including  the  free  surface, 
we  apply  successively  lower  order  of  the  central,  staggered  finite  difference  operators  used 
for  discretizing  the  spatial  derivatives.  OABC  is  used  at  the  bottom  and  sides  of  the  grids. 
It  is  worthwhile  to  note  that  in  the  discretization  we  had  to  define  the  variables  at  each  grid 
point  in  such  a  way  that  we  did  not  make  use  of  the  comers  of  the  computational  domain. 

ED  involves,  as  mentioned,  some  restrictions  on  the  computational  storage.  Applying 
this  method,  it  is  necessary  to  preserve  strips  along  the  boundaries  of  the  computational 
domain  and  to  introduce  damping  mechanisms  there.  The  exponential  damping  consists  of 
multiplying  each  physical  variable  in  the  strip  by  exp(-(ac  *  idist)2),  where  idist  is  the  cur¬ 
rent  number  of  grid  points  to  the  first  grid  point  of  absorption,  and  ac  is  a  constant  equal  to 
0.015  in  our  case.  Damping  is  enforced  within  absorption  strips  of  20  grid  points  along  the 


bottom  and  the  sides  of  the  grid,  while  gradual  decrease  of  the  spatial  finite  difference 
order  is  applied  only  when  approaching  the  free  surface.  In  both  methods,  boundary  con¬ 
ditions  as  for  a  plane  free  surface  are  used  at  the  upper  boundary  of  the  domain. 


Results 

In  our  numerical  tests  the  model  size  is  100  km  by  70  km,  with  351  grid  points  vertically 
and  501  grid  points  horizontally.  The  grid  length  in  both  directions  is  200  m.  A  pressure 
wave  is  initiated  from  a  Gaussian  point  source  with  a  central  frequency  of  2.5  Hz.  We 
chose  an  upper  bound  of  1.5%  for  the  relative  error  of  the  numerical  group  velocity.  This 
requires  3  samples  per  shortest  wavelength  in  case  of  the  8th  order  operator  used.  Further¬ 
more,  two  different  medium  models  are  tested;  a  homogeneous  one  where  the  P-velocity 
is  constant,  and  a  stratified  one  where  the  medium  P-velocities  are  constant  within  each 
layer. 

In  connection  with  our  use  of  higher  order  finite  differences  for  spatial  discretization,  we 
found  that  the  OABC  became  unstable  for  long  time  simulation  runs.  This  result  was  con¬ 
firmed  by  extensive  numerical  testing,  including  long  computer  runs,  i.e.  18376  time 
steps,  which  means  about  75  seconds  of  wavefield  synthetics.  Initially  we  used  8th  order 
of  accuracy  in  the  interior  of  the  domain  and  then  second  order  along  three  levels  next  to 
the  absorbing  boundaries  (thus  a  very  sharp  discontinuity  in  the  order  of  accuracy).  We  got 
instabilities  after  0(100)  time  steps;  instability  started  either  along  the  absorbing  bound¬ 
aries  or  at  the  corners  at  the  bottom  of  the  computational  grid.  From  there  it  spread  out  and 
destroyed  the  synthetics. 

Peng  and  Tokspz  (1994)  gave  some  examples  of  3D  elastic  finite  difference  runs  where 
OABC  was  used.  The  employed  scheme  had  fourth  order  accuracy  in  space  and  second 
order  accuracy  in  time.  Elastic  wave  propagation  in  an  unbounded  homogeneous  medium 
were  simulated,  and  the  results  were  stable. 
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Thus  for  a  realistic  test  of  OABC  using  our  scheme  it  was  natural  to  lower  the  spatial 
finite  difference  order.  That  choice  led  to  successful  runs  whether  using  2nd  order  finite 
difference  operators  all  over  the  grid  or  using  4th  order;  in  both  cases  we  got  completely 
stable  runs,  but  the  penalty  was  a  larger  numerical  dispersion. 

So  how  to  get  a  better  performance  for  our  higher  order  scheme,  in  implementing 
OABC,  in  terms  of  both  accuracy  and  stability?  We  solved  the  problem  by  decreasing 
gradually  the  finite  difference  order  when  approaching  the  boundaries.  It  was  a  matter  of 
experimentation  to  find  the  proper  number  of  grid  points  near  the  boundaries  on  which  to 
apply  a  certain  finite  difference  order.  It  seems  important  to  use  exactly  two  2nd  order  grid 
points  near  the  boundaries,  and  then  to  increase  the  number  of  grid  points  at  which  higher 
order  finite  differences  are  used. 

In  our  discretization  we  arrived  at  a  total  number  of  12  grid  points  for  which  to  lower 
gradually  the  finite  difference  order.  We  used  8th  order  in  the  interior  of  the  grid,  then  6th, 
4th  and  2nd  order  as  approaching  the  boundaries.  In  this  way  we  managed  to  run  for  at 
least  14000  time  steps  (about  60  seconds  of  simulated  time)  before  a  slowly  growing 
instability  started  appearing  along  the  absorbing  boundaries.  Simulated  wavefields  for 
such  large  time  laps  implies  that  the  most  interesting  wave  phenomena  are  adequately  syn¬ 
thesized. 

Peng  and  Tokspz  (1994)  mentioned  the  importance  of  stability  in  the  choice  of  an 
absorbing  boundary  condition.  As  an  example,  they  mentioned  that  Emerman  and  Stephen 
(1983)  demonstrated  that  the  boundary  condition  of  Clayton  and  Engquist  (1977)  was 
unstable  for  a  wide  range  of  elastic  parameters,  and  further  that  Mahrer  (1986)  found  the 
boundary  condition  of  Reynolds  (1978)  to  be  unstable  as  well.  Peng  and  Tokspz  made  a 
comparison,  using  a  4th  order  spatial  finite  difference  scheme  in  3D,  between  the  OABC, 
Reynolds’  and  Higdon’s  (1990)  boundary  conditions.  They  found  that  the  OABC  had  less 
artificial  reflections  than  the  two  latter  schemes.  Thus  our  study  may  be  considered  a  fur¬ 
ther  test  of  absorbing  boundary  conditions. 


In  order  to  compare  the  artificial  reflections  using  both  the  OABC  and  the  ED,  two  dif¬ 
ferent  types  of  models  were  used.  The  first  model  is  homogeneous  with  a  constant  P- 
velocity  of  7.1  km/s.  The  second  model  is  a  multilayer  realization  of  the  crust  and  upper 
mantle  consisting  of  several  constant  velocity  layers  (Fig.  5).  In  all  cases  the  S-wave  layer 
velocity  is  related  to  the  P-wave  velocity  via  Poisson’s  ratio;  by  this  choice  the  densities 
are  linear  functions  of  the  P-velocities.  At  the  top  of  the  model  free  surface  boundary  con¬ 
ditions  are  implemented.  A  Gaussian  point  source  is  located  50  km  from  the  left  edge  of 
the  grid  and  500  m  below  the  surface. 

For  the  first  two  runs  the  homogeneous  medium  model  was  used,  the  first  one  using  the 
OABC  (Figs  1  and  2)  and  the  second  one  the  ED  (Figs  3and  4).  In  Fig.  1  the  velocity  vec¬ 
tor  modulus  is  displayed  in  a  snapshot  11s  after  a  pressure  wave  has  been  initiated.  We  see 
that  the  P  wave  front  at  this  time  has  just  hit  the  bottom  of  the  model.  The  snapshot  reveals 
a  weak  P-P  reflection  from  the  bottom,  while  from  both  sidewalls  stronger  artificial  P-P 
and  P-S  reflections  can  be  noticed.  Interestingly,  the  reflections  from  the  left  sidewall  are 
stronger  than  the  ones  from  the  right  sidewall.  The  reason  for  this  is,  as  mentioned  by 
Levander  (1988),  that  we  make  use  of  two  staggered  numerical  space  differentiators  and  a 
staggered  grid  definition  of  the  field  variables.  The  consequence  is  that  not  exactly  the 
same  variables  are  being  absorbed  at  each  boundary.  Experiments  aimed  at  rectifying  this 
deficiency  of  asymmetric  reflections  were  done  by  using  the  OABC  to  absorb  every  vari¬ 
able  at  its  outermost  defined  point,  but  that  produced  even  stronger  grid  boundary  reflec¬ 
tions. 

In  Fig.  2  the  velocity  vector  modulus  is  displayed  in  a  snapshot  22  s  after  the  pressure 
wave  has  been  initiated.  We  see  mainly  artificial  reflections  here,  since  the  dominant  non¬ 
artificial  wavefield  has  passed  out  of  the  model  frame.  The  reflections  from  the  S-wave  are 
stronger  than  from  the  P-wave,  and  correspond  to  our  choice  of  maximum  P-wave  absorp¬ 
tion  (Peng  and  Toks0z,  1994).  The  strongest  reflections  are  seen  to  come  from  the  bottom, 
while  the  vertically  oriented  wave  front  seen  near  the  right  boundary  is  the  P-wave  reflec¬ 
tion  from  the  left  boundary. 


Figs  3  and  4  are  similar  to  Figs  1  and  2  respectively,  except  that  in  this  case  ED  was 
applied  within  strips  of  20  grid  points  along  the  numerical  boundaries  at  the  bottom  and 
sides.  Clearly,  the  ED  scheme  works  better  than  the  OABC  in  absorbing  both  P-waves  and 
S-waves.  It  absorbs  equally  well  along  the  bottom  and  the  sidewalls.  In  Fig.  4  the  non-arti¬ 
ficial  waves  have  passed  out  of  the  model,  and  only  artificial  grid  reflections  remain.  They 
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are  weaker  than  reflections  with  the  OABC  by  a  factor  of  approximately  10  ,  and  fur¬ 
thermore,  the  run  is  completely  stable. 

A  homogeneous  model  is  rather  an  oversimplification  of  the  real  Earth,  and  so  further 
comparison  between  OABC  and  ED  is  tied  to  wavefield  simulation  using  a  multilayered 
crust  -  upper  mantle  model  shown  in  Fig.  5.  Velocity  vector  modulii  as  in  Figs  1-4  are 
shown  in  Figs  6-9,  but  now  with  the  background  model  of  Fig.  5  added.  This  fact,  together 
with  the  multiply  additional  reflections  from  the  various  layer  boundaries,  cause  these 
snapshots  to  be  less  clear  than  those  for  the  homogeneous  medium  model  (Figs  1-4).  Figs 
6  and  7  show  the  wave  field  respectively  11  seconds  and  22  seconds  after  the  pressure 
wave  was  initiated,  using  OABC.  Also  in  this  case  we  see  asymmetric  P-P  and  P-S  reflec¬ 
tions  from  the  sidewalls,  together  with  a  whole  range  of  wave  conversions  and  reflections 
from  the  various  layers.  The  absorption  at  the  bottom  is  seen  to  be  somewhat  better  than 
that  at  the  sidewalls.  This  is  due  to  the  staggered  grid  discretization. 

Figs  8  and  9  show  corresponding  runs  where  we  made  use  of  ED.  We  see  the  same  wave 
conversions  and  reflections  as  those  of  Figs  6  and  7  caused  by  the  layered  medium,  but 
now  less  prominent.  Along  the  numerical  boundaries  the  absorption  strips  of  20  grid 
points’  width  are  clearly  visible.  We  get  a  confirmation  again  that  the  performance  of  ED 
is  better  than  that  of  OABC.  Comparing  the  results  of  the  last  two  runs  we  found  that  the 
artificial  reflections  from  using  OABC  were  on  the  average  0(10)  stronger  than  the  ones 
using  ED.  Noteworthy,  for  the  long  runs  using  ED,  the  damping  affects  the  resolution  in 
the  interior  of  the  domain,  slightly  weakening  the  wave  amplitudes.  In  this  respect  OABC 


behaves  better. 


Conclusions 


In  implementing  the  OABC  we  found  the  choice  of  spatial  finite  difference  order  to  be 
important.  It  must  be  decreased  when  approaching  the  numerical  boundaries  in  order  to 
keep  the  run  stable  as  long  as  possible.  Our  synthetic  discretization  scheme  uses  8th  order 
finite  differences  in  the  interior  of  the  grid.  Thus,  moving  from  the  interior  towards  the 
boundaries  we  found  the  following  procedure  to  be  optimal:  6  layers  of  6th  order  finite 
differences,  then  4  layers  of  4th  order  and  finally  2  layers  of  2nd  order  next  to  the  numeri¬ 
cal  boundary.  Comparing  the  amount  of  numerical  boundary  reflections  using  OABC  in 
relation  to  using  ED,  it  seems  clear  that  ED  works  better,  even  when  applying  it  within  a 
thin  strip.  The  mentioned  drawback  of  additional  computer  storage  are  offset  by  the 
smaller  simulation  time  using  ED. 
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FIGURE  1.  Snapshot  of  velocity  vector  modulus  11  s  after  a  pressure  wave  is  initiated  from  a 
Gaussian  point  source  at  depth  500  m  below  the  surface.  The  medium  model  is  homogeneous  and 
OABC  is  used. 
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FIGURE  2.  Snapshot  of  velocity  vector  modulus  for  the  same  situation  as  in  Fig.  1, 11  seconds  later. 
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FIGURE  3.  Snapshot  of  velocity  vector  modulus  for  the  same  situation  as  in  Fig.l,  but  in  this  case 
ED  is  used. 
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FIGURE  4.  Snapshot  of  velocity  vector  modulus  for  the  same  situation  as  in  Fig.  3, 11  seconds  later. 
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FIGURE  5.  Snapshot  of  velocity  vector  modulus  13  s  after  a  pressure  wave  is  initiated  from  a 
Gaussian  point  source  at  depth  500  m  below  the  surface.  The  source  position  and  the  background 
model  of  the  layered  medium  are  clearly  visible.  The  maximum  P-velocity  is  8.18  km/s  right  below 
the  Moho.  The  Moho  is  located  at  35  km  depth.  The  nearest  layer  above  the  Moho  represents  a  P- 
velocity  of  6.9  km/s.  The  next  layer  is  the  maximum  velocity  layer  in  the  crust  with  a  P-velocity  of 
7.1  km/s.  Just  below  the  surface  there  is  the  minimum  velocity  layer  with  a  P-velocity  value  of  6.6 
km/s. 
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FIGURE  6.  Snapshot  of  velocity  vector  modulus  11  s  after  a  pressure  wave  is  initiated  from  a 
Gaussian  point  source  at  depth  500  m  below  the  surface.  The  medium  model  is  layered  and  OABC 
is  used. 
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FIGURE  7.  Snapshot  of  velocity  vector  modulus  for  the  same  situation  as  in  Fig.  6, 11  seconds  later. 
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FIGURE  8.  Snapshot  of  velocity  vector  modulus  for  the  same  situation  as  in  Fig.6,  but  in  this  case 
ED  is  used. 
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FIGURE  9.  Snapshot  of  velocity  vector  modulus  for  the  same  situation  as  in  Fig.  8, 11  seconds  later. 
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1.  Introduction 

The  transformation  of  seismic  profiling  data  to  the  domain  of  intercept  time 
and  ray  parameter  (r  —  p)  and  the  subsequent  interpretation  and  inversion 
for  velocity-depth  models  by  the  r-sum  method  have  become  a  popular 
processing  technique  in  wide-angle  marine  seismic  profiling  (Stoffa  et  al., 
1981;  Diebold  and  Stoffa,  1981).  The  method  offers  several  advantages: 

—  In  processing: 

The  r  -  p  transformation  is  a  beamforming  process  (plane  wave  de¬ 
composition)  in  which  coherent  events  are  enhanced  and  incoherent 
ones  suppressed. 

—  In  interpretation: 

The  r-sum  inversion  method  require  that  we  identify  and  digitize  the 
‘principal  branch’  in  the  r-p  domain  consisting  mainly  of  post-critical 
reflections  and  diving  waves.  Since  post-critical  reflections  are  totally 
reflected  they  have  relatively  large  amplitudes  compared  to  first  ar¬ 
rivals  which  usually  are  head  waves.  Thus  the  problem  of  picking  weak 
first  arrivals  (which  is  often  a  critical  process  when  interpreting  and 
inverting  data  in  x  -  t  domain)  is  avoided. 

—  In  Inversion: 

The  simple  r-sum  traveltime  inversion  technique  of  Diebold  and  Stoffa 
(1981)  is  easy  to  implement  and  works  both  for  discontinuous  and 
continuous  velocity-depth  functions.  Its  main  limitation  is  that  the 
velocity  must  increase  monotonically  with  depth. 
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Like  many  other  inversion  methods  r-sum  is  limited  to  laterally  homoge¬ 
neous  models,  i.e.,  an  one-dimensional  model  in  the  sense  that  the  velocity 
is  a  function  of  depth  only.  The  case  of  dipping  plane  layer  interfaces  is 
discussed  by  Diebold  and  Stoffa  (1981).  In  this  paper  we  will  examine  the 
effect  of  relatively  small  scale  inhomogeneities  introduced  in  the  form  of 
random  perturbations  with  a  von  Karman  autocorrelation  function  and 
an  RMS  of  3%.  We  generate  synthetic  data  by  a  2-dimensional  8-th  order 
elasto-dynamic  finite  difference  method  (Hestholm  et  ah,  1994;  Sguazzero 
et  ah,  1989)  for  both  a  laterally  homogeneous  and  a  perturbed  model.  The 
synthetic  seismograms  are  transformed  to  r  —  p  domain,  interpreted,  and 
velocity- depth  functions  are  computed  for  both  models.  Besides  comparing 
the  two  resulting  models  we  can  also  compare  the  data  in  both  x  —  t  and 
r  -  p  domain. 


2.  Synthetic  model  and  seismogram  generation 

The  P-wave  velocity-depth  function  of  the  unperturbed,  laterally  homoge¬ 
neous  model  is  given  in  Fig.  1.  The  crustal  model  consists  of  both  constant 
velocity  and  gradient  layer.  The  S-wave  velocities  were  calculated  by  as¬ 
suming  a  constant  Poison’s  ratio  of  0.25.  Densities  were  found  from  Birch’s 
law. 

For  the  finite  difference  computation  we  used  a  340  km  long  and  70  km 
deep  model.  The  grid  step  is  0.2  km  and  the  source  signature  is  a  Ricker 
wavelet  with  a  center  frequency  of  2.5  Hz  (the  data  were  bandpass  filtered 
1  to  4.5  Hz  prior  to  processing).  Absorbing  boundaries  are  used  at  the 
sides  and  bottom  of  the  model,  while  vanishing  stresses  are  used  at  the  free 
surface.  The  explosion  source  is  placed  0.5  km  below  the  surface  10  km 
from  the  left  side  of  the  model.  The  seismograms  of  60  seconds  duration 
are  extracted  at  800  points  along  the  surface  to  the  right  of  the  source  and 
with  0.4  km  intervals. 

The  perturbed  model  was  generated  by  scaling  the  unperturbed  one 
with  a  2D  random  field  having  a  unit  mean,  a  von  Karman  autocorrela¬ 
tion  function  of  order  0.3  (corresponding  to  a  fractal  dimension  of  2.7), 
correlation  distances  of  10  km  horizontally  and  2.5  km  vertically,  and  an 
RMS  value  of  0.03  (3%).  The  generation  and  use  of  stochastic  field  in  seis¬ 
mic  modelling  is  discussed  in  e.g.,  Frankel  and  Clayton  (1981),  Charrette 
(1991),  and  Holliger  and  Levander  (1992).  The  choice  of  these  model  param¬ 
eter  is  discussed  in  Hestholm  et  al.  (1994),  and  is  believed  to  be  realistic  in 
a  crystalline  continental  crust  environment.  A  section  from  our  perturbed 
model  is  shown  in  Fig.  2. 
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3.  Finite  difference  modelling  results 


The  synthetic  seismogram  profiles  are  shown  in  Fig.  3  and  4  for  the  unper¬ 
turbed  and  perturbed  models  respectively.  In  both  cases  the  dominating 
phase  (highest  amplitude)  is  the  fundamental  mode  Rayleigh  wave  (Rg) 
propagating  with  a  group  velocity  of  3.15  km/s.  For  small  offsets  (say  less 
than  30  km)  the  seismogram  section  from  the  unperturbed  model  is  rel¬ 
atively  easy  to  interpret.  It  consists  mainly  of  the  direct  P-wave,  the  Rg 
phase  and  P  reflections  from  the  model  discontinuities.  With  increasing  off¬ 
set  we  see  a  more  complicated  P  wavetrain  (coda)  consisting  of  P-multiples 
and  mixed  waves,  i.e.,  waves  that  have  propagated  partly  as  P  and  partly  as 
S  due  to  mode  conversion  at  internal  reflectors  or  at  the  free  surface.  These 
mixed  waves  will  have  apparent  velocities  in  the  same  range  as  pure  P-waves 
because  in  a  laterally  homogeneous  model  the  ray  parameter  (slowness)  will 
not  change  during  reflections  or  transmission.  At  distances  beyond  about 
145  km  the  Moho  refracted  Pn  phase  is  the  first  arrival,  but  of  relatively 
weak  amplitude. 

The  ideal  explosion  source  used  in  this  modelling  experiment  does  not 
generate  S-waves  when  applied  in  a  homogeneous  unbounded  medium,  but 
because  the  source  here  is  close  to  the  surface  (compared  to  wavelength) 
quite  strong  S-waves  are  generated  in  all  directions  except  for  the  vertical. 
The  S-coda  looks  simpler  than  the  P-coda  in  this  profile  because  in  the  dis¬ 
tance  range  out  to  about  150  km  (where  we  have  a  opportunity  to  observe 
the  first  part  of  the  S-coda)  the  multiple  S-phases  have  too  high  appar¬ 
ent  velocities  to  be  trapped  in  the  crustal  waveguide.  However,  at  larger 
distances  the  S-multiples  (constituting  the  Lg  phase)  have  lower  velocities 
(<  4.5  km/s)  and  thus  propagate  as  guided  waves. 

The  seismic  profile  from  the  perturbed  model  (Fig.  4)  is  very  messy 
compared  to  Fig.  3.  The  P-wave  reflections  from  the  discontinuities  in  the 
original  model  are  difficult  to  identify  and  looks  like  coda  waves  even  at 
short  offsets.  The  Pn  phase  shows  a  variable  amplitude,  it  is  stronger  than 
for  the  unperturbed  model  in  some  distance  ranges  and  almost  disappear 
at  others.  The  P-coda  is  strong  into  the  S-waves  and  masks  most  of  the  Sn 
phase.  Much  of  the  coda  waves  are  seen  to  consist  of  waves  propagating 
with  Rg  velocities  and  in  both  forward  and  backward  directions.  These 
scattered  surface  waves  are  generated  when  a  strong  primary  wave  hits 
strong  heterogeneities  near  the  free  surface.  Also  in  real  recordings  surface 
waves  are  seen  to  constitute  a  significant  part  of  the  coda  (see  e.g.,  Hestholm 
et  al.,  1994). 
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4.  r  —  p  transformation 

The  subsequent  processing  of  the  synthetic  dataset  was  performed  using 
the  ProMAX  commercial  processing  system.  Prior  to  r  —  p  transformation 
a  f  —  k  filter  was  applied  in  order  to  remove  the  dominating  Rg-wave.  The 
fan-shaped  f  —  k  filter  rejected  all  waves  with  apparent  velocities  less  than 
3.33  km/s.  With  the  small  intersensor  interval  of  0.4  km  for  this  profile, 
spatial  aliasing  is  not  a  problem  for  the  frequencies  of  interest,  i.e.,  less  than 
4.5  Hz.  In  order  to  avoid  artifacts  in  the  transformed  data  due  to  truncation 
effects  we  also  found  it  necessary  to  mute  the  data  both  in  time  and  space. 
A  5  sec  long  ramp  function  was  used  at  the  end  of  each  seismogram,  and 
similarly  a  20  km  long  ramp  function  was  used  to  gradually  weight  down 
the  seismograms  at  the  end  of  the  profile. 

The  resulting  r  —  p  transformed  datasets  are  shown  in  Fig.  5  and  6 
for  the  unperturbed  and  perturbed  model  respectively.  As  for  the  time- 
offset  profiles  there  is  a  great  difference  between  the  models.  Although 
Fig.  5  is  rather  complicated  with  both  P-,  S-  and  mixed- waves,  the  different 
phases  forms  distinct  branches  in  r  —  p  domain  and  most  of  these  are  easily 
explained  after  some  consideration.  As  expected  the  reflected  P-waves  forms 
quarter-ellipses  in  the  upper  left  hand  corner  of  the  figure  (r  range  0  -  11  s 
and  p  range  0  -  0.165  s/km).  The  corresponding  S-wave  branches  are  found 
by  multiplying  both  the  r  and  p  values  by  the  P/S-velocity  ratio  of  1.73.  The 
S- waves  are  very  clear  in  the  p  range  0.165  to  0.285  s/km  where  only  S- waves 
can  exist  in  this  model.  For  slownesses  below  0.1  s/km  the  S- waves  gradually 
disappear  because  the  source  does  not  generate  vertically  propagating  S- 
waves.  The  most  complex  part  of  the  figure  is  for  the  slowness  range  0.125 
to  0.165  s/km  for  which  the  crust  acts  as  a  partial  waveguide  for  P-waves 
(some  P-wave  energy  will  leak  through  the  Moho  when  converted  to  S- 
waves  at  the  free  surface  or  at  model  discontinuities).  The  ray  parameters 
bounding  this  interval  is  also  seen  to  be  slownesses  at  which  abrupt  changes 
in  the  waveform  occur  (due  to  phase  shifts  in  the  reflection  coefficients). 
The  lower  value  (0.125  =  1/8.0)  is  determined  by  the  P-velocity  below 
the  Moho  and  the  upper  value  (0.165  ^  1/6.1)  by  the  P-velocity  near  the 
surface.  A  similar  change  in  waveforms  is  seen  in  the  Moho  reflected  S-wave 
at  the  ray  parameter  corresponding  to  the  S-velocity  below  the  Moho.  The 
rather  strong  wave  with  an  intercept  time  of  about  12  sec  and  a  slowness  of 
about  0.10  s/km  is  the  Moho  P-to-S  reflection.  Its  amplitude  depends  on 
the  reflection  coefficient  which  have  a  maximum  for  slownesses  around  0.08 
s/km.  Similar  P-to-S  reflections  are  visible  also  for  some  internal  reflectors 
in  the  crust.  In  addition  to  the  phases  discussed  above  some  multiple  P- 
and  S-waves  are  seen  below  the  'pure’  P-  and  S-branches.  In  the  slowness 
interval  0.125  to  0.165  many  of  the  phases  are  probably  multiples  of  mixed 
(P-S)  waves. 
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The  r  -  p  transformed  data  for  the  perturbed  model  (Fig.  6)  is  quite 
messy  as  compared  to  Fig.  5  and  not  unexpected  when  comparing  the  two 
datasets  in  the  time-offset  domain.  The  part  that  is  most  easily  recognized 
when  compared  to  Fig.  5  is  the  post-critical  P-reflections.  Except  from 
this  only  some  parts  of  the  pre-critical  S-reflections  and  possibly  the  Moho 
reflected  P-phase  is  reasonably  intact.  All  other  phases  are  wiped  out  by 
interfering  and  scattered  wavelets. 


5.  r-sum  inversion 

Fortunately,  the  post-critical  P-reflections  that  are  most  suitable  to  invert 
for  the  velocity  depth  model  is  also  the  part  of  the  r  —  p  data  that  is  best 
preserved  when  introducing  model  perturbations.  The  digitized  parts  of 
the  ‘principal  branch’  is  shown  in  Fig.  7  for  both  models.  The  results  are 
surprisingly  consistent  considering  the  otherwise  remarkable  differences  in 
the  two  data  sets.  There  might  be  some  ambiguity  in  choosing  which  cycle 
to  sample  but  in  this  case  the  same  cycle  appears  to  have  been  chosen  for 
both  models.  The  resolution  is  limited  by  the  signal  bandwidth.  In  real 
wide-angle  reflection/refraction  profiling  work  signals  are  recorded  up  to 
much  higher  frequencies  (say  10-20  Hz)  than  those  used  in  our  synthetic 
experiments.  However,  with  the  rather  high  degree  of  scattering  usually 
observed  in  continental  crystalline  crust  at  high  frequencies,  the  effective 
frequency  range  (limited  by  the  spatial  correlation)  is  not  much  different 
from  that  used  here. 

The  digitized  r  —  p  curves  were  inverted  to  velocity-depth  functions  by 
means  of  the  r  —  sum  method  (Diebold  and  Stoffa,  1981).  Here,  a  constant 
velocity  layer  is  assigned  to  every  sample  of  the  T-p  curve.  By  dense  r  -p 
sampling  one  can  achieve  an  almost  continuous  velocity- depth  function.  It  is 
not  necessary  to  know  if  any  part  of  the  curve  represent  reflections  or  diving 
waves  in  a  gradient  layer.  If  there  is  some  ‘noise  in  the  sampled  curve, 
artifacts  like  negative  layer  thicknesses  may  occur,  but  will  not  severely 
corrupt  the  end  results.  The  inherent  assumption  of  no  low  velocity  zones 
in  the  model  is  no  problem  in  our  experiment.  The  resulting  velocity- depth 
function  for  both  models  and  the  original  input  model  is  shown  in  Fig.  8. 
The  general  trend  of  the  model  is  well  recovered  for  both  the  perturbed 
and  unperturbed  dataset.  None  of  the  two  recovered  models  is  significantly 
better  than  the  other.  In  both  cases  the  differences  between  the  obtained 
model  and  the  original  is  less  than  0.3  km/s  for  almost  all  depths.  For  both 
models  the  obtained  velocities  are  generally  too  low  in  the  upper  part  of  the 
crust  (down  to  about  20  km).  This  might  be  a  coincidence,  but  probably 
reflect  difficulties  in  sampling  the  r  -  p  curve  at  small  intercept  times  (see 
Figs.  5  and  6). 
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6.  Conclusions 

Although  most  features  of  the  r  -  p  transformed  data  for  the  unperturbed 
model  was  lost  for  the  perturbed  model,  we  were  able  to  use  the  post- 
critical  reflections  in  a  r  —  sum  inversion  to  recover  for  both  models,  with 
about  the  same  accuracy,  the  essential  features  of  the  original  velocity-depth 
distribution.  Futher,  for  both  models  the  finer  details  of  the  layering  is  lost 
and  the  resolution  of  the  obtained  velocity-depth  functions  seems  to  be 
about  0.2  km/s.  The  resolution  is  limited  by  the  usable  signal  bandwidth 
which  again  depend  on  the  degree  of  scattering.  In  our  experiment  the 
dominating  frequency  of  the  signal  was  about  2.5  Hz  and  the  RMS  of  the 
perturbed  medium  was  3%  for  a  von  Karman  correlation  function  of  order 
0.3.  Velocity-depth  models  obtained  by  r  —  p  transformation  and  r  —  sum 
inversion  seems  to  give  a  good  laterally  averaged  velocity  in  the  presence 
of  velocity  heterogeneities  of  the  size  used  in  this  experiment,  i.e.,  about  10 
km,  and  recording  offsets  up  to  about  300  km. 
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Figure  1.  The  P-velocity  of  the  unperturbed  model.  The  Moho  is  at  31  km.  The  same 
velocity-depth  function  was  used  as  a  basis  for  the  perturbed  model  (Fig.  2). 
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Figure  2.  A  section  showing  the  P-velocities  of  the  perturbed  model.  The  (horizon¬ 
tally)  averaged  velocities  are  the  same  as  for  the  unperturbed  model  (Fig.  1).  The  2D 
auto-correlation  function  of  the  velocity  perturbation  field  is  a  von  Karman  function  of 
order  0.3,  horizontal  and  vertical  correlation  distances  of  10  km  and  2.5  km  respectively, 
and  an  RMS  of  3%. 
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Figure  3.  The  time-offset  section  generated  from  the  unperturbed  model.  Maximum 
offset  (horizontal  axis)  is  320  km.  Time  (vertical  axis)  in  seconds.  The  amplitude  scaling 
is  constant  for  the  whole  section. 
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Figure  4 .  The  time-offset  section  generated  from  the  perturbed  model.  Maximum  offset 
(horizontal  axis)  is  320  km.  Time  (vertical  axis)  in  seconds.  The  amplitude  scaling  is 
constant  for  the  whole  section. 


10 


B.O.  Ruud  ET  AL. 


Figure  5.  The  r  —  p  transformed  data  from  the  unperturbed  model.  Slowness  p  along 
horizontal  axis  in  s/km.  Intercept  time  r  along  vertical  axis  in  seconds.  See  text  for 
explanation  of  the  different  branches. 
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Figure  6.  The  t  -  p  transformed  data  from  the  perturbed  model.  Slowness  p  along 
horizontal  axis  in  s/km.  Intercept  time  r  along  vertical  axis  in  seconds. 
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Figure  7.  The  ‘principal  branch5  of  the  r  —  p  transformed  data  digitized  from  Figs.  5 
and  6. 


55  6  6.5  7  7.5  8  8.5 

Velocity  (km/s) 


Figure  8.  The  velocity-depth  functions  obtained  by  r-sum  inversion  of  the  r  —  p  curves 
shown  in  Fig.  7.  For  comparison  the  original  model  is  also  included. 
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The  general  popularity  of  magnitude  as  a  convenient  and  robust  measure  of  earthquake  size  makes  it 
tempting  to  examine  whether  this  parameter  can  be  reliably  estimated  in  near  real  time.  In  this  study 
we  demonstrate  that  this  is  indeed  the  case  conditioned  on  the  design  of  the  signal  detector  being  of 
STA/LTA  type  where  STA  is  a  short  term  signal  power  or  rms  estimate.  Using  real  data  we  demon¬ 
strate  the  Random  Vibration  Theory  relation  that  /!max  -  (2\nN)^  A  m5  is  valid  for  non-stat.onary 
seismic  signals.  Using  Rayleigh’s  theorem  we  also  established  a  relation  between  X™,  and  the  flat  por¬ 
tion  of  the  source  spectra.  These  Amax  and  Arms  estimation  procedures  are  used  for  determining  con¬ 
ventional  magnitudes  and  moment  magnitudes  for  29  events  as  recorded  by  the  Norwegian  Seismo¬ 
graph  Network  (NSN).  We  used  here  a  procedure  outlined  by  Sereno  et  a/.  (1988)  and  also  their  geo¬ 
metrical  spreading  and  attenuation  parameters  derived  from  analysis  of  NORSAR  recordings.  Our 
magnitude  and  moment  magnitude  estimates  for  5  different  frequency  bands  are  m  good  agreement 
with  the  Ml  estimates  derived  from  the  conventional  magnitude  formulas  in  combination  with  empin 
cal  correction  tables.  Surprisingly,  the  Hmax  and  magnitudes  produced  consistent  negative  bmsed 
bv  ca.  0.4  units  estimates  even  in  the  extreme  4-8  Hz  band.  In  view  of  the  good  agreement  between 
various  tvpes  of  magnitude  estimates,  we  constructed  conventional  magnitude  correction  tables 
spreading  and  attenuation  parameters  from  Sereno  et  al  (1988)  for  a  variety  of  signal  frequency 
bands.  Near  real  time  Amax  ad/or  Arms  or  correspondingly  event  magnitudes  would  be :  of ^ significance  m 
automatic  phase  association  analysis,  bulletin  production  for  local  and  regional  seismic  networks  and 
the  earthquakes  monitoring  performances  of  such  networks. 


Key  words  random  vibration  theory  -  distance 
corrections  —  spreading  and  attenuation  effects  - 
network  monitoring  -  capabilities  in  real  time 


1.  Introduction 

The  concept  of  earthquake  magnitude 
was  firstly  introduced  by  Richter  (1935). 
He  proposed  a  logarithmic  amplitude  scale 
tied  to  the  first  arriving  P  wave  as  a  mea¬ 
sure  for  the  relative  size  of  local  Southern 
California  earthquakes.  This  simple  and  ro¬ 
bust  magnitude  scale  being  popular  with 


seismologists  around  the  world,  it  has  been 
extended  to  all  distance  ranges  and  many 
seismic  phases  (Nuttli,  1973;  Bath  et  al., 
1981;  Ebel,  1982).  Despite  the  diversity  of 
record  measurements,  tied  to  the  maximum 
amplitudes  of  P,  S  and  surface  waves,  all 
such  scales  are  of  similar  forms 

M  =  log  +  q(rJ)  +  c  (I-*) 

where  M  is  the  magnitude,  A  the  maximum 
phase  amplitude,  T  is  the  period  (s),  q  ac- 
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counts  for  geometrical  spreading  and  atten¬ 
uation  as  a  function  of  distance  r  (km)  and 
signal  frequency  (f)  while  c  is  a  station  cor¬ 
rection  term.  Focal  depth  (h)  is  seldom  in¬ 
cluded  explicitly  in  magnitude  formulas  al¬ 
though  it  is  significant  in  case  of  surface 
waves  from  deep  events. 

Although  an  occurring  earthquake  is  a 
complex  phenomenon,  the  Richter  magni¬ 
tude  remains  highly  useful  in  many  seismo- 
logical  contexts.  A  drawback  with  such 
measurement  using  conventional  instru¬ 
mentation  is  that  earthquake  phase  ampli¬ 
tudes  may  be  measured  at  frequencies  be¬ 
yond  the  corner  frequency  of  the  source 
spectra  which  cause  a  consistent  magnitude 
bias.  Therefore,  many  seismologists  argue 
in  favour  of  replacing  magnitude  with  Aki’s 
(1967)  earthquake  moment  M0  which  may 
be  obtained  from  inversion  of  broadband 
waveform  data  (Dziewonski  and  Gilbert, 
1974;  Jost  and  Herrmann,  1989)  or  equiva¬ 
lently  the  level  of  the  flat  portion  of  the 
source  spectra.  To  avoid  confusion,  mo¬ 
ment-magnitude  scales  have  been  intro¬ 
duced  in  the  form  of  linear  relationships 
between  the  moment  (M0)  and  magnitude. 
For  example,  Sereno  et  al.  (1988)  intro¬ 
duced  the  following  formula  for  local  Scan¬ 
dinavian  earthquakes: 

log  (M0)  =  1.03A4  +  17.1  (1.2) 

where  ML  is  the  local  magnitude. 

In  eq.  (1.1)  the  AIT  measurement  re¬ 
mains  essentially  a  manual  operation  while 
the  g-term  is  derived  from  empirical  tables 
when  an  epicenter  solution  is  at  the  hand. 
In  this  paper  we  address  the  problem  of  un¬ 
dertaking  in  near  real  time  event  magni¬ 
tude  and/or  moment  measurements.  This 
study  is  tied  to  analysis  of  local  events  be¬ 
cause  A/L-magnitude  estimation  is  often 
problematic  due  to  complex  geometrical 
spreading  and  attenuation  relationships  at 
small  distance  ranges.  It  is  considered  an 
important  task  in  view  of  the  widespread 
use  of  event  magnitude  in  observational 
seismology  and  the  potential  of  using  am¬ 
plitude  information  in  automatic  schemes 


for  epicenter  location  (Ruud  et  al ,  1993; 
Sambridge  and  Gallagher,  1993)  and  net¬ 
work  performance  simulations  (Sereno  et 
al ,  1988). 


2.  On  line  AJT  estimation 

Obviously  for  near  real  time  magnitude/ 
moment  estimation,  the  A/T  term  must  be 
an  integral  part  of  the  signal  detector  de¬ 
sign.  The  most  popular  ones  in  use  are  the 
so-called  sliding  window  STA/LTA  type 
which  is  a  comparison  of  short  and  long 
term  trace  amplitude  (a,)  averages  (e.g., 
Ruud  and  Husebye,  1992).  Common  STA 
definitions  are  of  the  forms: 


STA  (abs)  =—2  |  at  |  (2.1) 

n  i  =  l 


STA  (rms) 


1/2 


(2.2) 


Equivalent  expressions  are  valid  for 
LTA  when  the  effective  window  length  n  is 
increased  10-20  times  that  of  the  STA  win¬ 
dow  length.  The  STA  and  LTA  terms  are 
commonly  measured  in  near  real  time  on 
bandpass  filtered  seismometer  output 
traces.  When  the  STA/LTA  ratio  exceeds  a 
present  threshold,  the  presence  of  a  seismic 
signal  is  declared.  Our  preference  is  for  the 
rms  related  STA-LTA  definitions  also  be¬ 
cause  signal  rms  is  related  to  the  signal 
power  spectrum  via  the  Parsevaks  theo¬ 
rem. 

From  eq.  (2.2)  the  task  is  to  estimate 
peak  time  domain  amplitudes  for  the  first 
arrival  of  P- waves  or  the  dominant  Lg - 
phase  in  case  of  local  events  that  is  the 
A  IT  At  rm  in  eq.  (1.1).  From  the  Random 
Vibration  Theory  (RVT)  results  of 
Cartwright  and  Longuet-Higgins  (1956)  we 
have  the  relation 


E{Am3X)  =  ArmJ{N)  (2.3) 


where  E  is  expectation,  N  is  the  number  of 
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extremes  (peaks  and  troughs)  and  >lnns  is 
the  rms  value  of  the  signal.  In  the  case  of 
STA  trace  window,  Arms  -  STA(rms).  We 
have  for  /(TV) 

/(TV)  =  (21n  TV)172  (2.4) 

The  RVT  assumes  stationary  time  series 
(Crandal  and  Mark,  1963)  which  are  sel¬ 
dom  the  case  for  seismic  recordings.  How¬ 
ever,  Boore  (1983)  and  Boore  and  Joyner 
(1984)  have  demonstrated  that  the  RVT  re¬ 
lation  in  eq.  (2.3)  performs  well  even  in 
cases  of  peak  signal  accelerations  and  ve¬ 
locities  as  determined  from  time  domain 
simulations.  The  validity  of  the  simple  lin¬ 
ear  relationship  in  eq.  (2.3)  was  confirmed 
by  extensive  analysis  of  real  data,  and  such 
results  are  detailed  in  the  data  analysis  sec¬ 
tion. 

Whether  we  prefer  to  measure  magni¬ 
tude  (Mi)  or  moment  (M0),  the  Amax  term 
is  estimated  from  bandpass  filtered  traces, 
and  thus  must  be  corrected  for  the  instru¬ 
ment  response.  A  reasonable  assumption 
here  is  that  the  displacement  response  spec¬ 
tra  are  approximately  flat  within  the  filter 
passband,  and  that  the  instrument  response 
correction  is  tied  to  the  center  frequency  of 
the  bandpass  filter.  Similarly,  the  period 
(T)  of  the  maximum  signal  amplitude  is 
also  tied  to  the  center  frequency  of  the 
bandpass  filter.  Thus,  both  the  Amax  and  T 
terms  in  magnitude  formula  (eq.  (1.1))  are 
obtained  via  eq.  (2.3). 


3.  Magnitude  measurements 

The  AIT- term  derived  via  the  detector 
parameter  STA(rms)  is  sufficient  for  event 
magnitude  estimation  given  that  a  tabula¬ 
tion  of  the  distance  dependent  g-term  in 
eq.  (1.1)  is  available.  However,  for  many 
regions  the  g-term  is  not  well  known  at  lo¬ 
cal  and  regional  distances  making  corre¬ 
sponding  event  magnitude  estimates  unreli¬ 
able.  For  example,  as  a  rule  the  ISC  does 
not  estimate  P-wave  magnitudes  for  r  < 
20°.  In  other  cases,  available  correction  ta¬ 


bles  were  derived  from  old  analog  record¬ 
ings  in  the  1-2  Hz  bandpass  while  modern 
instrument  recordings  of  local  events  are 
dominated  by  3-6  Hz  signal  frequencies. 
An  additional  problem  is  that  intraplate 
seismicity  levels  are  so  low  that  it  is  difficult 
to  calibrate  the  local  event  magnitudes  to¬ 
wards  those  obtained  from  teleseismic 
recordings. 

In  other  words,  local  magnitude  scales 
may  differ  significantly  from  each  other,  so 
recently  seismologists  have  been  exploring 
the  usefulness  of  event  moment  magnitude 
scales. 


3.1.  Moment  measurements 


As  mentioned  above,  the  q-term  in  eq. 
(1.1)  is  not  always  well  established,  and  a 
way  of  circumventing  this  problem  is  to  in¬ 
troduce  specific  estimates  of  the  geometri¬ 
cal  spreading  and  attenuation  ( Q )  effects. 
Following  Sereno  et  al.  (1988),  we  write  for 
the  amplitude  spectrum  of  displacement  at 
epicenter  distance  r 


|  A {f,  r)  |  =  S{f)G{r,  r0)  exp 


-71  ft  ' 

eco. 


(3.1) 


where  /  is  frequency,  S(f)  source  spectra, 
G(r,  r0)  geometrical  spreading  function 
with  a  reference  distance  (r0),  and  Q(f )  the 
attenuation  function.  The  attenuation  func¬ 
tion  is  of  the  form 

Q(f)  =  QoF  (3-2) 

where  the  frequency  dependence  is  through 
the  fy  expression.  The  geometrical  spread¬ 
ing  term  is  of  the  form 

f  (1  /rY  ;  r  <  rQ 

G(r,  r0)  =  (3-3) 

l  (Vro)(ro/ry  ;  r  >  r0 

where  the  rQ,  ft,  and  v  parameters  depend 
on  the  phase  type  and  besides  may  vary 
from  one  tectonic  regime  to  another.  At  lo- 
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cal  distances  the  dominant  signal  in  the 
records  is  the  Lg-phase  whose  source  spec¬ 
trum  in  case  of  an  earthquake  is  expressed 
as  (Sereno  et  al ,  1989) 

S(f)  -  SqH(J,  fc)  (3.4) 


H(f,  /c)  = 

=(1  +  (1  -*‘>(iY+*{£YY'B 

(3.5) 


The  source  spectrum  decays  with/-2  be¬ 
yond  the/c  due  to  the  H  function.  From  eq. 

(3.6)  we  can  estimate  the  SQ  an d/c  parame¬ 
ters.  However,  in  practice  the  corner  fre¬ 
quency  is  not  well  resolved,  so  our  prefer¬ 
ence  was  to  fix  it  at  a  «safe»  value  of  20  Hz 
which  is  also  very  close  to  the  bandwidth  of 
the  signals  used  in  analysis  (sampling  fre¬ 
quency  =  50  Hz).  Under  this  assumption, 
the  H(f fc)  is  negligible  and  combining 
eqs.  (3.1),  (3.6),  (3.7)  and  (3.8)  we  get 

log  M0  (max)  ^ 


where  50  is  the  low  frequency  spectral 
source  level,  fc  is  the  corner  frequency  of 
the  source  spectrum,  B  is  the  amount  of 
overshoot  (Xie,  1993).  We  have  the  follow¬ 
ing  relations  for  fc  and  S0  regarding  M0: 


fc  — 


1/3. 


So 


M0 

Anofi3 

(3.6) 


-  log  An 


log  2A fi  +  log  4rtop3  + 


log  G(r,  r0)  +  ■ 


Tift 

W) 


log  e  (3.9) 


log  Mo(rms)  “ 

=  log^rms  +  ~T  log  -^~r  +  log  4jto(33 
z  ~a.fi 


where  /3  is  the  crustal  shear  wave  velocity,  g 
the  crustal  density  and  c  a  scaling  constant. 
For  moment/magnitude  estimation,  we 
must  establish  a  relationship  between  the 
amplitude  spectrum  of  the  displacement 
(|  A(f,r)  |  )  and  the  near  «real-time»  record 
parameters,  or  Amzx,  namely 

I  A  (foh  r)  I  =  fA^rms  ;  i  =  1,...,5 
r  2A fi 

(3.7) 

I  ^  (foi,  f)  \  —  0  Af  ^  max  >  ^  1,...,5 

(3.8) 

where  fQi  represents  the  center  frequencies 
of  the  5  bandpass  filters  used  in  data  analy¬ 
sis,  |  A(f0i,  r)  |  is  the  amplitude  spectrum  at 
frequency  f0i  and  distance  r,  AT  is  the  win¬ 
dow  length  in  s  and  A f  is  the  bandwidth  of 
the  rlh  filter.  The  derivations  of  eq.  (3.7) 
and  eq.  (3.8)  are  detailed  in  Appendix. 


-  log  G(r,  r0)  +  ~fL.\oge  (3.10) 

Note  that  P9  S  and  Lg-phases  from  ex¬ 
plosions  and  earthquakes  respectively  have 
different  source  representations.  However, 
the  main  point  is  that  it  should  be  feasible 
to  estimate  in  near  real  time  event  magni¬ 
tude  and  event  moment  from  «real-time» 
available  signal  detector  parameters.  The 
validity  of  this  statement  will  be  demon¬ 
strated  in  the  next  section. 


4.  Data  analysis  and  results 

In  this  section  we  would  demonstrate  that 
the  deduced  Amax  -  A ^  relationship  in  eq. 
(2.3)  appears  to  be  valid  for  the  recordings  of 
local  earthquakes  and  explosions  at  stations 
equipped  with  short  period  seismometers. 
We  use  observations  from  the  Norwegian 
Seismograph  Network  (NSN)  whose  station 
locations  are  shown  in  fig.  1.  A  more  subtle 
but  extensive  test  was  performed  by  calculat- 
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Fig.  1.  Stations  constituting  the  Norwegian  Seismograph  Network  (NSN)  which  is  operated  by  the 
University  of  Bergen.  Some  of  the  stations  are  equipped  with  three-component  instrumentation.  In 
our  analysis  only  the  vertical  short  period  recordings  are  used.  The  NORSAJR  operated  regional  ar¬ 
rays  NORESS  (NAO,  NRAO)  and  ARCESS  (ARAO)  are  also  marked. 
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Table  I.  Listing  of  all  the  events  used  in  our  analysis.  The  focal  parameters  are  taken  from  the  mon¬ 
thly  bulletins  for  the  Norwegian  Seismic  Network  (NSN)  as  published  by  the  University  of  Bergen. 
The  Mc  and  ML  notations  reflect  NSN  coda  wave  magnitude  and  NORSAR  reported  ML  (L^-waves) 
magnitude,  respectively.  We  have  introduced  a  few  modifications  in  parameter  listing;  all  presumed 
explosions  (EXP)  are  given  zero  focal  depth.  The  separation  of  events  into  EXP  and  EQ  (earthquake) 
populations  are  not  based  on  specific  classification  criteria.  Potential  errors  here  would  not  significan¬ 
tly  change  our  analysis  results. 


N. 

Date 

Time 

Location 

Mc 

H 

(km) 

# 

Stat. 

Source 

type 

(d/m/y) 

(h.min:s) 

(LatN) 

(LonE) 

1 

26/08/93 

10.34:32.0 

59.05 

5.77 

2.6 

2.2 

0 

8 

EXP 

2 

26/08/93 

19.22:04.9 

61.11 

4.07 

2.8 

2.5 

18 

8 

EQ 

3 

10/09/93 

13.11:23.8 

58.34 

6.36 

2.7 

2.0 

0 

5 

EXP 

4 

13/09/93 

05.25:14.4 

66.34 

5.67 

3.9 

3.2 

11 

10 

EQ 

5 

15/09/93 

09.58:38.1 

60.60 

4.79 

2.1 

1.6 

0 

7 

EXP 

6 

15/09/93 

15.31:29.7 

67.11 

20.84 

3.2 

1.8 

0 

4 

EXP 

7 

21/09/93 

13.15:21.1 

58.33 

6.31 

3.0 

2.5 

0 

7 

EXP 

8 

28/09/93 

16.44:53.9 

58.49 

10.62 

2.9 

2.3 

15 

5 

EQ 

9 

28/09/93 

16.57:42.1 

60.46 

5.15 

1.6 

1.4 

0 

6 

EXP 

10 

29/09/93 

14.12:21.5 

58.15 

6.31 

2.8 

2.2 

0 

6 

EXP 

11 

02/10/93 

08.30:55.3 

60.39 

5.02 

1.8 

1.5 

0 

6 

EXP 

12 

03/10/93 

23.48:35.5 

60.06 

7.29 

2.5 

1.8 

15 

7 

EQ 

13 

04/10/93 

20.21:48.3 

61.94 

1.44 

2.3 

1.0 

11 

5 

EQ 

14 

07/11/93 

23.40:44.0 

67.84 

20.08 

2.5 

1.7 

15 

4 

EQ 

15 

07/11/93 

23.43:17.2 

66.28 

7.02 

3.1 

2.5 

25 

4 

EQ 

16 

08/11/93 

14.28:17.3 

59.91 

2.72 

2.4 

1.5 

15 

7 

EQ 

17 

12/11/93 

19.54:36.0 

59.69 

12.86 

3.4 

2.9 

16 

6 

EQ 

18 

15/11/93 

15.04:25.6 

62.17 

3.27 

2.7 

1.9 

15 

5 

EQ 

19 

21/11/93 

01.53:56.3 

60.18 

4.96 

2.6 

2.1 

12 

6 

EQ 

20 

27/11/93 

18.57:52.7 

60.47 

11.66 

3.1 

2.2 

0 

5 

EXP 

21 

13/12/93 

09.00:09.5 

56.75 

2.74 

3.1 

2*6 

12 

6 

EQ 

22 

27/12/93 

05.20:46.3 

61.29 

2.79 

3.3 

3.6 

14 

8 

EQ 

23 

03/01/94 

22.12:59.8 

61.76 

4.19 

2.6 

2.4 

17 

7 

EQ 

24 

07/01/94 

09.06:22.6 

60.60 

2.43 

2.3 

1.8 

23 

5 

HQ 

25 

15/01/94 

00.00:26.2 

65.23 

7.58 

2.8 

1.6 

2 

3 

EQ 

26 

19/01/94 

09.16:54.8 

66.36 

14.60 

2.4 

1.4 

0 

3 

EXP 

27 

21/01/94 

01.26:14.1 

65.98 

11.89 

3.1 

2.1 

19 

2 

EQ 

28 

25/01/94 

03.07:58.3 

62.46 

5.07 

2.6 

1.9 

20 

5 

EQ 

29 

26/01/94 

17.27:47.1 

66.84 

13.58 

3.0 

2.3 

0 

5 

EXP 

ing  magnitudes  and  moments  for  the  29 
events  listed  in  table  I  using  the  approach  de¬ 
tailed  in  the  previous  section. 


4.1.  Amax  deduced  from  ATms 

In  figs.  2a-f  to  4a-f  the  apparent  validity 
of  the  Amax~ATms  relationship  is  demon¬ 
strated  for  station  KMY  recording  of  Event 
4  in  table  I.  In  (a)  the  P  and  Lg  phase 


recordings  are  shown;  (b)-(c)  the  indepen¬ 
dently  observed  Amax  versus  Arms  or  STA 
for  both  P  and  Lg  waves;  (d)-(e)  the  corre¬ 
sponding  number  of  extremes  within  the  in¬ 
dividual  STA  windows.  Note  that  an  ex- 
trem  occurs  when  the  second  derivative  of 
the  trace  time  function  is  zero.  For  unfil¬ 
tered  traces,  as  in  figs,  (d)  and  (e),  the 
number  of  extreme  would  be  relatively  high 
due  to  the  ripples  overlaying  the  signal 
trace.  The  linear  Amax  -A^  relationship  in 
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e)  Lg  Phase 
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Fig.  2a-f.  Amax-Arms  analysis  for  KMY  unfiltered  recordings  of  Event  4  in  table  II  at  an  epicentral  di¬ 
stance  of  795  km.  a)  Original  vertical  trace  record;  b)  and  c)  Am„  versus  Arm$  for  P-  and  Z.?-phases 
respectively.  The  window  length  used  in  the  evaluation  of  ^4rms  is  2.5  s;  d)  and  e)  number  of  extremes 
iV  for  P-  and  L^-phases;  f)  observed  and  estimated  Amax  using  eq.  (2.3).  The  two  curves  match  each 
other  almost  perfectly. 
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Filter :  2.0  -  4.0  Hz 
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Fig.  4a-f.  Amax-Ams  analysis  for  KMY  recordings  in  the  2.0-4.0  Hz  passband  of  Event  4  in  table  II 
Otherwise,  caption  as  in  fig.  2a-f. 
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figs,  (b)  and  (c)  can  be  predicted  from  the 
number  of  extremes  (N)  in  figs,  (d)  and 
(e);  the  agreement  is  very  good.  From  this 
result  we  expect  an  excellent  agreement  be¬ 
tween  observed  and  estimated  Amax  values 
which  indeed  is  the  case  from  results  pre¬ 
sented  in  figs.  (f).  Even  in  the  extreme 
cases  of  very  peaked  Z^-waves,  this  agree¬ 
ment  remains  good. 

So  far,  we  have  computed  the  N- param¬ 
eter  for  every  STA  window:  since  this  pa¬ 
rameter  fluctuates  moderately  within  a  spe¬ 
cific  filter  band  and  the  (21n  N)1/2  expres¬ 
sion  of  eq.  (2.3)  vary  slowly  with  N ,  there 
is  really  no  need  to  update  the  //-estimates 
for  every  STA  window.  Anyway,  the  essence 
of  the  results  presented  in  figs.  2a-f  to  4a-f  is 
that  the  STA  values  as  estimated  in  near  real¬ 
time  for  a  rms  type  of  STA/LTA  detector  can 
be  used  for  accurately  estimating  the  maxi¬ 
mum  trace  amplitudes  for  any  part  of  the 
short  period  seismometer  records. 


4.2.  Magnitude  estimation 

Given  the  epicentral  distance  and  Amax 
or  ^rmsj  we  can  compute  event  magnitudes 
or  event  moments  as  detailed  above.  This 
was  done  for  the  29  events  listed  in  table  I 
and  then  these  moments  were  converted 
to  event  magnitudes  using  the  conversion 
formulas  and  the  parameter  values  given 
in  table  III.  These  moment  magnitudes 
are  denoted  A/L(Ser)max  and  A/L(Ser)rmr 
Since  frequency  dependence  on  A/L-esti- 
mates  is  a  debated  topic,  we  computed 
ML{ Ser)  for  5  frequency  bands  as  indicated 
in  table  II.  Using  the  maximum  amplitude 
values  for  S  or  L^-phases  we  also  computed 
conventional  event  magnitudes  using  vari¬ 
ous  types  of  correction  curves  (spreading 
and  attenuation)  as  taken  from  NORSAR, 
Alsaker  et  al.  (1991),  and  the  Seismological 
Observatory  in  Helsinki  (details  also  in 
table  III).  In  the  latter  cases,  there  is 
hardly  any  frequency  dependence  in  the 
distance  correction  term,  so  here  only  the 
f2  (1.5-3  Hz)  frequency  band  was  used  in 
the  magnitude  calculations. 


In  table  II  all  magnitude  values  repre¬ 
sent  averages  with  respect  to  the  number  of 
reporting  stations  listed  in  table  I.  The  ob¬ 
tained  event  magnitudes  are  similar  at  the 
lower  frequencies,  irrespective  of  epicenter 
distance.  The  only  exception  here  is  the 
MLh  which  are  consistently  high  by  ca.  0.6 
magnitude  units.  A  comparison  with  the 
coda  magnitudes  (Mc)  gives  that  these  are 
also  consistently  high  and  besides  appear  to 
be  less  stable  than  the  NORSAR  reference 
magnitudes  (ML).  The  latter  are  in  reason¬ 
able  agreement  with  the  /2  ML(Ser)  values 
in  table  II  that  are  consistently  low  by  ca. 
0.3  magnitude  units.  These  features  are  il¬ 
lustrated  in  fig.  5a, b  showing  plots  of 
A/L(Ser)max  versus  ML-magnitudes  re¬ 
ported  by  NORSAR. 

Frequency  dependence  of  the  magnitude 
estimates  is  obvious  from  the  table  II  re¬ 
sults,  but  there  is  hardly  any  difference  be¬ 
tween  the  ML(Ser)rms  and  A/L(Ser)max  esti¬ 
mates.  The  attenuation  parameter  may  be 
important  in  this  respect,  so  we  recalculated 
the  ML( Ser)  magnitudes  for  y  =  0.20  and 
y  ~  0.14  equivalent  to  increased  signal  at¬ 
tenuation.  Its  effect  is  to  increase^the  /4 
and  /5  magnitudes  by  only  0.1  magnitude 
units  in  average  since  epicenter  distances 
are  generally  less  than  250  km.  On  the 
other  hand,  the  differences  between /I  and 
/2  as  compared  to  /4  and  f5  magnitude  es¬ 
timates  remain  at  ca.  0.4  magnitude  units. 
In  fig.  6a-f  we  have  plotted  event  magni¬ 
tudes  for  individual  stations  as  a  function  of 
distance  for  y  =  0.26,  y  =  0.40  and  y  = 
0.55  which  are  equivalent  to  decreasing  at¬ 
tenuation.  In  these  plots  the  average  event 
magnitude  is  adjusted  to  a  value  of  2.5 
units,  while  the  stations  at  distances  less 
than  250  km  are  deleted  since  the  effect  of 
frequency  dependence  is  as  mentioned  very' 
small  for  such  ranges.  From  fig.  6a-f,  we 
see  that  a  y  value  of  0.40  or  0.55  would  be  a 
better  choice  than  y  =  0.26.  These  results 
contradicting  those  in  table  II.  are  not  con¬ 
sidered  conclusive  since  we  have  few  obser¬ 
vations  beyond  500  km  and  the  scattering 
in  the  individual  stations  ML  values  is  of 
the  order  of  some  tens  of  a  magnitude  unit. 
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Table  II  Magnitude  estimates  for  the  events  listed  in  table  I.  fl, ft, ft, ft  and/5  represent  dlff®rel^ 
lame  11.  „  36  ^  and  4_g  Hz>  respectively.  Column  ML  is 

L  6  e  ve  n  t^m a°ni  tu de s  as  reported  by  NORSAR  using  NORESS  recordings  while  Mc  is  coda  magnitude 
based  NSN  recordings  as  reported  in  the  bulletin.  A/L(Ser)max  is  magnitude  estimates  tied  to  maxi¬ 
mum  amplitude  observations  (eq.  (3.8))  while  A/L(Ser)tms  are  corresponding  moment  magnitude  esti¬ 
mates  tied  to  the  rms  values  (eq  (3.7)).  Common  features  are  averaging 

n  and  usine  Sereno  et  al  (1988)  spreading  and  attenuation  constants  (table  III).  Columns  MLN,  MLA 
aid  ml  '  am magnitude  estimates  simila°r  to  ML(Ser)max  except  that  the  distance  correction  terms 
a  nf  Alsaker  et  al  (1991)  and  Helsinki  (table  III).  Since  these  are  frequency 

ta  fheMld  „e  !Ld,  Note  .he  <imU.riqr  of  even,  m.gntodes 

inlhe  ft  band  with  the  exception  of  ML H  and  Mc  estimates.  Another  notable  feature  is  the  consistent 
decay  in  ML(Ser)  from  the  fl-ft  bands  towards  the /4-/5  bands.  We  take  this  to  imply  that  the  SILg 
excitation  is  less  effective  at  higher  signal  frequencies. 


N.  Ml 


ML( Ser)max 

~fl  ft  ft  ~ft 


ft  ft 


ML  (Ser)rms 

~ft  ft  ~ft 


MLn  MLa  MLh 

~ft~  ~ft  ft  ~ft 


2.6  2.5 

2.4  2.3 


The  above  attenuation  results  signify  the 
importance  of  path  effects  and  station 
corrections,  and  also  the  observational 
fact  that  Lg- wave  excitation  is  relatively  less 
effective  at  signal  frequencies  above  3-4 
Hz. 


5.  Discussion 

In  this  study  we  have  explored  the  possi¬ 
bility  of  estimating  event  magnitude  in  near 
real*  time.  At  this  stage  of  development  we 
have  concentrated  on  demonstrating  that 
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Table  III.  Listing  of  correction  formulas  as  used  for  the  estimated  magnitudes  tabulated  in  table  II 
M0  =  moment,  A  =  peak  amplitude  in  nm  and  A  =  epicentral  distance  in  km.  F( A,  T)  correction 
function  is  shown  in  fig.  7  for  T  =  0.45  s. 


Magnitude  formulations 

Reference 

log  M0  =  1.03  ML{ Ser)  -  17.1  (EQ) 
log  M0  =  1.04  ML  (Ser)  -  17.7  (EXP) 

Qo  =  560;  y  =  0.26 
r0  =  100  km;  a  =  1;  v  =  1/2 

Q  =  2.7  gr/cm3;  =  3.5  km/s 

Sereno  et  ai  (1989) 

MLn  =  log  (100  A)  +  F( A,  T) 

NORSAR 

MLa  =  log  A  +  0.91  log  A  +  0.00087  A  -  1.31 

Alsaker  et  ai  (1991) 

MLh  =  log  A  +  1.27  log  A  -  1.44 

Helsinki  Seismological  Obser/atory 

Fig.  sa.b.  ML(Ser)max  magnitude  estimates  from  table  II  versus  the  reference  NORSAR  ML  magni¬ 
tude  estimates  in  table  II.  The  small  offset  may  in  part  reflect  local  site  effects  since  the  NORESS 
array  is  located  outside  the  NSN  siting  area. 
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Amax  and/or  ^4^  and  hence  event  magni¬ 
tudes  can  be  obtained  from  signal  detector 
outputs  of  the  STA/LTA  types.  In  figs.  2a-f 
to  4a-f,  an  excellent  agreement  between 
observed  and  estimated  maximum  ampli¬ 
tude  in  all  frequency  bands  considered  is 
demonstrated.  However,  we  have  not  con¬ 
templated  such  maximum  amplitude  esti¬ 
mates  in  a  real  time  operational  environ¬ 
ment  since  we  are  not  directly  involved  in 
this  kind  of  work. 

Through  the  analysis  of  29  events  (earth¬ 
quakes  and  explosions)  we  have  deter¬ 
mined  y4max  and^4rms  as  described  above  for 
both  conventional  and  moment  magnitude 
estimates.  In  the  former  cases,  we  used  em¬ 
pirical  formulas  for  the  geometrical  spread¬ 
ing  and  attenuation  effects  taken  from 
Alsaker  et  al.  (1991)  and  those  used  at 
NORSAR  and  Helsinki  in  their  bulletin 
works.  In  many  regions  such  correction  ta¬ 
bles  are  not  available,  so  we  explored  the 
possibility  of  using  moment  magnitudes 
based  on  spreading  and  attenuation  param¬ 
eters  estimated  directly  from  spectral  analy¬ 
sis  of  local  event  records.  All  the  various 
types  of  Ml  estimates  for  local  events 
recorded  by  the  Norwegian  Seismograph 
Network  (NSN)  stations  proved  to  be  mu¬ 
tually  consistent.  The  only  exception  here 
was  the  ML  Helsinki  estimates  which  ap¬ 
peared  to  be  positive  biased  by  ca.  0.6  mag¬ 
nitude  units.  In  the  NSN  monthly  bulletins, 
coda  magnitudes  are  routinely  reported. 
These  estimates  also  appear  to  be  positive 
biased  and  somewhat  unstable. 

Initially,  we  allowed  the  corner  fre¬ 
quency  (fc)  to  be  a  free  parameter  but  for 
larger  events  ML  >  3.0,  this  caused  erro¬ 
neous  estimates  in  the  sense  that  unrealistic 
magnitudes  estimates  were  obtained.  The 
remedy  here  was  to  freeze  the  corner  fre¬ 
quency  at  20  Hz.  This  also  implies  that  the 
estimate  of  the  fc  parameter  is  not  well  con¬ 
strained  in  the  analysis  procedure  used.  For 
larger  events,  say  for  ML  >  4.5,  some  of 
the  frequency  bands  considered  here  would 
be  beyond  the  comer  frequency  and  hence 
the  corresponding  moment  magnitude  esti¬ 
mates  should  have  negative  biases.  Indica¬ 


tions  of  such  cases  would  be  negative  dif¬ 
ferences  between  low  and  high  frequency 
bands,  actually  observed  for  nearly  all 
events  used  in  analysis  (table  II).  From  the 
previous  section,  we  rule  out  relatively 
strong  signal  attenuation  or  exceptionally 
low  corner  frequencies  as  a  plausible  expla¬ 
nation  for  these  observational  features.  In¬ 
stead,  we  consider  that  the  S/Lg  signal  exci¬ 
tation  spectra  are  not  flat  met  to  the  corner 
frequency  as  presumed  in  our  theoretical 
modelling.  In  other  words,  the  most  consis¬ 
tent  and  correct  magnitude  estimates  are 
obtained  for  signal  frequencies  in  1-3  Hz 
range;  at  higher  signal  frequencies  the  ob¬ 
served  negative  biases  can  be  computed  by 
introducing  frequency  dependent  station 
corrections  if  considered  worth  while. 

A  recent  tendency  in  establishing  local 
magnitude  scales  is  the  synthesis  of 
the  original  Wood-Anderson  seismograph 
recordings  from  Southern  California 
(Richter,  1935)  (e.g.  Alsaker  et  al. ,  1991). 
Our  preference  is  for  the  moment  magni¬ 
tude  estimation  approach  simply  because  of 
its  anchoring  on  physical  source  representa¬ 
tion  and  well  established  wave  propagation 
parameters.  These  parameters  can  be  ob¬ 
tained  through  analysis  of  appropriate  sig¬ 
nal  spectra  (Sereno  et  al ,  1988)  which  be¬ 
sides  would  be  representative  of  the  sta¬ 
tion/network  siting  area.  In  contrast,  the 
wave  propagation  regime  of  Southern  Cali¬ 
fornia  is  not  representative  for  most  of  the 
network  areas  elsewhere  and  hence  not 
well  suited  as  a  magnitude  reference  base. 
In  a  seismic  hazard  study  for  Southern  Nor¬ 
way,  Singh  et  al.  (1990)  discussed  various 
magnitude  formulas  in  use  and  also  the  as¬ 
sociated  attenuation  parameters  including 
those  of  Sereno  et  al.  (1988).  The  differ¬ 
ences  here  appear  larger  than  the  corre¬ 
sponding  differences  in  estimated  event 
magnitudes  and  a  reconciling  explanation 
here  is  that  spreading  and  attenuation  ef¬ 
fects  are  not  easily  separated.  The  com¬ 
bined  effects  are  incorporated  in  the  magni¬ 
tude  correction  curves.  The  relative  consis¬ 
tency  of  AZL(Ser)  magnitudes  (table  II)  im¬ 
plies  that  it  should  be  possible  to  construct 
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Fig.  6a-f.  The  ML(Ser)  estimates  as  a  function  of  epicenter  distance  and  different  values  of  the  atte¬ 
nuation  parameter  y.  The  corresponding  linear  regression  lines  are  also  shown.  Note  that  network 
A/L(Ser)max  average  for  each  event  is  scaled  to  a  reference  magnitude  of  2.5  units  and  this  scaling  con¬ 
stant  would  be  different  for  each  frequency  band  considered.  Observations  for  distances  below  250  km 
were  excluded  from  the  figure.  There  are  very  few  observations  beyond  500  km  so  it  is  not  quite  ob¬ 
vious  that  the  best  fitting  y  =  0.40  (diminishing  attenuation)  is  truly  correct.  The  table  II  results  give 
less  differences  between  low  and  high  frequency  estimates  for  the  decreasing  y  values. 
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Fig.  7.  Generating  distance-correction  curves  from  the  Sereno  et  al  (1988)  spreading  and  attenuation 
parameters  as  listed  in  table  III.  For  comparison  the  corresponding  Alsaker  et  al  (1991)  and  NORSAR 
(7  =  0.45  s)  curve  are  shown.  The  correction  curves  derived  from  Sereno  et  al  (1988)  are  plotted  for 
the  filters /l,/3  and/5.  In  table  II  the  ML(Ser)  and  MLA  magnitude  estimates  are  very  similar  which 
is  also  expected  from  fig.  7  since  epicenter  distances  generally  are  less  than  250  km.  The  NORSAR 
reported  magnitudes  (ML  column  in  table  II)  are  in  average  0.3  units  lower  than  expected  since  the 
corresponding  epicenter  distances  are  in  the  300-600  km  range.  Note  that  the  different  correction  cur¬ 
ves  are  similar  out  to  ca.  350  km  that  is  the  distance  range  where  the  majority  of  local  events  are 
recorded. 


correction  tables,  the  q(r,  /)-term  in  eq. 
(1.1),  from  the  spreading  and  attenuation 
parameters  in  eq.  (3.9).  Such  correction 
curves  are  shown  in  fig.  7  for  the  center  fre¬ 
quencies  and  bandwidths  of  the  bandpass 
filters  used  and  the  Adsaker  et  al.  (1991) 
and  NORSAR  correction  curves  are  in¬ 
cluded  for  comparison.  In  table  II  the  MLA 
magnitudes  are  very  similar  to  the  ML(Ser) 
ones  as  expected  from  fig.  7  since  the  epi¬ 


center  distances  seldom  exceed  250  km. 
Likewise,  the  NORSAR  reported  NtL  mag¬ 
nitudes  are  in  average  0.3  units  below  the 
ML(Ser)  ones.  This  is  also  easily  explained 
from  fig.  7  since  the  NORSAR  epicenter 
distances  (see  fig.  1)  are  mostly  in  the  300- 
600  km  range  where  the  associated  correc¬ 
tion  curve  is  lower  by  ca.  0.3  units.  In  other 
words,  our  magnitude  estimates  and  dis¬ 
tance  corrections  pertain  to  epicenter  dis- 
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tances  well  below  500  km.  Ruud  and  Huse¬ 
bye  (1992)  have  demonstrated  that  at  least 
in  Fennoscandia  most  of  the  station  event 
recordings  pertain  to  epicenter  distances 
less  than  500  km.  If  magnitude  is  to  be  re¬ 
tained  as  a  standard  earthquake  size  mea¬ 
sure,  a  universal  moment  to  magnitude  for¬ 
mula  may  be  constituted  say  by  IASPEI,  in 
order  to  avoid  the  prevailing  confusion  on 
type  of  event  magnitude  as  published  for 
local  and  regional  events.  A  final  remark  is 
that  near  real  time  event  magnitude  esti¬ 
mates  would  be  invaluable  for  estimating 
network  monitoring  capabilities  also  in 
near  real  time.  Such  estimates  are  tied  to 
noise  level  estimates  from  the  individual 
network  stations  in  combination  with  a  grid 
search  procedure  ( e.g .  Sereno  and  Bratt, 
1989;  Ringdal  and  Kvaerna,  1992).  Another 
application  now  under  consideration  is  to 
use  near  real  time  magnitude  estimates  as  a 
part  of  automated  epicenter  location 
schemes  (Ruud  et  al. ,  1993;  Sambridge  and 
Gallanger,  1993)  that  is  by  incorporating 
both  phase  and  amplitude  consistency  in 
the  phase  association  routine  and  hence  in 
the  final  epicenter  location. 


6.  Concluding  remarks 

In  this  study  we  have  demonstrated  that 
near  real  time  magnitude  estimation  is  fea¬ 
sible  conditioned  on  a  signal  detector  de¬ 
sign  for  determining  peak  signal  trace  am¬ 
plitudes.  The  results  presented  in  table  II 
imply  that  moment  magnitudes  are  a  viable 
alternative  to  conventional  magnitude  esti¬ 
mates.  The  advantages  are  an  anchoring  to 
seismic  source  theory  and  the  use  of  clearly 
defined  signal  spreading  and  attenuation 
parameters.  We  have  also  demonstrated 
that  the  spreading  and  the  attenuation  pa¬ 
rameters  can  be  used  for  construction  of 
correction  tables  commonly  used  in  con¬ 
ventional  magnitude  determinations,  thus 
avoiding  tedious  analysis  of  many  earth¬ 
quake  recordings  (Alsaker  et  al. ,  1991). 
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Appendix  _ _ _ _  . 

The  relation  between  and  )A(f,  r)  I : 

Rayleigh’s  theorem  gives  the  energy  relation  between  the  time  and  frequency  domain  representations  such 
as: 

J”  |  a(t)  \2dt  =|”  \A(f)  \2df  (A.1) 

where  a(t)  is  the  displacement  in  time  and  \A(f)  |  is  the  amplitude  spectrum.  The  unbiased  rms  estimate  of 
the  signal  a(t)  in  continuous  time  is  defined  as: 

ALs  =-Tr-  J  I  a( 0  \2dt  (A.2) 

where  AT  is  window  length  in  s.  Likewise,  for  digital  signals  we  have 


N 


I  «? 


(A3) 


where  /V  is  the  number  of  the  samples.  From  eq.  (3.1),  we  assume  that  \  A  (J,  r)  |  is  flat  within  the  pass- 
band  of  the  filters  used  and  this  in  combination  with  the  Rayleigh’s  theorem  leads  to: 


A2  = 

■Ai  rms 


1 

AT 


1 2df  + 


A(f)  Vdf 


1 

AT 


2 A/  |  A(f0)\2 


(A.4) 
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where  fB  is  the  bandwidth  of  a(0>/i  and/2  are  the  low  and  high  cutoff  frequencies  of  the  bandpass  filter,  A f 
is  the  filter  bandwidth  (A f  =  f2  —  /i)  and/0  is  the  center  frequency  L/o  =  ifi  +  />)/2]  of  the  bandpass  filter. 
Thus,  eq.  (A.4)  gives  us  the  relationship  between  the  rms  value  of  the  signal  a(t)  and  the  amplitude  of  the 
spectrum  at  the  frequency  /0  which  is  expressed  as: 


I  A(foh  r ) 


(A.5) 


where /o/  and  A/  are  the  center  frequencies  and  the  bandwidths  of  the  5  bandpass  filters  used  in  data  analy¬ 
sis,  | Affoij  r)|  is  the  amplitude  spectrum  of  the  displacement  at  epicentral  distance  r.  In  data  analysis  AT 
being  the  STA  window  length  is  fixed  at  5  s  independent  of  the  bandpass  filter  used.  In  the  context  of  de¬ 
tector  design  the  STA  window  lengths  seldom  exceed  a  few  seconds,  however  it  is  easy  to  combine  short 
STA(rms)  values  in  order  to  produce  longer  rms. 


The  relation  between  Amax  and  jA(f,  r) j: 

Assuming  that  the  amplitude  spectrum  of  the  displacement  a{t)  is  flat  within  an  idealized  bandpass  filtered 
trace,  we  have: 

I  A(f)\  =  AQTI±f(f  -  /0)  +  A0TiAf(f  +  f0)  (A.6) 

where  nA/-  is  a  rectangular  function  with  unit  amplitude  and  bandwidth  A ftf0  is  the  center  frequency  of  nAy 
and  Aq  is  the  amplitude  of  the  displacement  spectrum.  If  A(f)  is  a  real  function,  its  inverse  Fourier  trans¬ 
form  is: 

a(t)  =  2AfA0sinc(Aft)cos(2jrfQt)  (A.7) 

a(t)  has  a  maximum  value  at  t  =  0,  namely 

Amax  =  2AfA0 

=  2Af  |  A(foi,  r)  |;  /  =  1  5  (A.8) 


Eq.  (A.7)  is  only  valid  in  case  of  real  functions  (no  phase  shifts)  which  is  hardly  the  case  in  seismology. 
However,  since  the  sine  function  varies  rather  slowly  as  compared  to  the  cos  function  in  eq.  (A.7),  it  can 
be  shown  that  the  Amax  relation  in  eq.  (A.8)  remains  valid  for  practical  purposes  also  for  non-real  func¬ 
tions. 


